CODEWITHSUNDEEP
r
Munish
Munish

Reputation: 677

Overlaying world map on a image plot in R

I am trying to plot a global map using latitude, longitude and grid data in R. For this I am using image and image.plot functions. Additionally I need to overlay global coastline for land area. However I am not sure how to place the map exactly over my image of gridded data. Map is appearing bit shifted to the left in the console and that part is not visible either. See sample code below with random grid data.

remove(list=ls())

library(fields)

library(maps)

grid_lon<-c(0.5:1:359.5)

grid_lat<-c(-89.5:89.5)

temp1<-matrix(data = rexp(200, rate = 10), nrow = 360, ncol = 180)#random matrix

zlim=c(0,0.25)

par(oma=c( 3,0,0,0))# c(bottom, left, top, right)#plot margins

image(grid_lon,grid_lat,temp1,axes=FALSE,xlab='',ylab='')

map("world", fill=TRUE, col="white", bg="white", ylim=c(-90, 90),add=TRUE)

title(main ='Main title')
image.plot(zlim=zlim,legend.only=TRUE,horizontal=TRUE,legend.mar=0.4,legend.shrink=0.4,legend.width=0.4,nlevel=64,axis.args = list(cex.axis =1,at=zlim, labels=zlim,mgp=c(1, 0, 0),tck=0),smallplot=c(.25,.72, 0,.030),
    legend.args=list( text=expression(textstyle(atop('anomaly',
    paste('(meters)')),cex.main=1.2)),cex=1.2, side=1, line=1.6)
    )#end image.plot

box()

Upvotes: 1

Views: 3825

Answers (2)

Munish
Munish

Reputation: 677

I found the answer after few attempts and a tip from colleague. What needs to be done is shift the longitude grid from 0:359 to -179.5:179.5 using following commands after grid_lon is declared: 

indexes_to_shift<-180

grid_lon[grid_lon>=180]<-grid_lon[grid_lon>=180]-360

grid_lon<-c(tail(grid_lon, indexes_to_shift), head(grid_lon, indexes_to_shift))

Upvotes: 0

Peter Herman
Peter Herman

Reputation: 162

In general, when working with maps it is preferable to use spatial objects, for which a projection method can be defined. The coherence with the map is then better guaranteed. Since you are working with a filled grid, an obvious choice is to use a raster from package raster. Your code would then become:

require (raster)
require (maps)
temp1<-matrix(data = rexp(180*360, rate = 10), nrow = 360, ncol = 180)  #random matrix
r<-raster(temp1,xmn=-179.5,xmx=179.5,ymn=-89.5,ymx=89.5,crs="+proj=longlat +datum=WGS84")
plot(r)
map("world",add=T,fill=TRUE, col="white", bg="white")

EDIT

This code does not take into account that the data comes as a 360*180 matrix, while it is desirable to plot (map) a 180*360 matrix. Transposing is risky because it may result in an upside-down image. In order to be sure that the right coordinates are associated with the right values, we can explicitly associate them, and afterwards transform into a spatial object. The for-loop doing this is in the code below is slow, maybe it can be made more efficient, but it does the job.

require (raster)
require (maps)
# basic data, as in code given
grid_lon<-seq(0.5,359.5,1)
grid_lat<-seq(-89.5,89.5,1)
temp1<-matrix(data = rexp(200, rate = 10), nrow = 360, ncol = 180)#random matrix
# transform into data frame, where coords are associated to values
tt<-data.frame(lon=rep(NA,64800),lat=rep(NA,64800),z=rep(NA,64800))
ct<-0
for (i in 1:360){
  for (j in 1:180){
    ct<-ct+1
    tt$lon[ct]<-grid_lon[i]
    tt$lat[ct]<-grid_lat[j]
    tt$z[ct]<-temp1[i,j]
  }
}
# transform to spatial structure
coordinates(tt)<- ~lon+lat
# make spatial structure gridded
gridded(tt)<-TRUE
# transform to raster
r<-raster(tt)
projection(r)<-crs("+proj=longlat +datum=WGS84")
# plot
plot(r)
map("world2",add=T,fill=TRUE, col="white", bg="white")

Upvotes: 1

Related Questions