Masking the desired country in spatial plot

[monamidutta]

Hi I am using esd to create field object and spatial plot which looks quiet ok. Since I have selected a rectangular domain for India (68,98,8,38), the data also includes some parts of neighbouring countries and the plot is given in attached file named plot1, however I want to mask the other countries (Pakistan, Tibet etc) and wanted a plot like plot2. I have tried map((mask(x),type="fill") which is masking only the ocean. Thanks in advance

[DIncongnito]

Hello there! That is a nice question. You can do it easily with a tiny trick. I hope you know the regrid() function :D. If you have an IMD gridded data, you can make a spatial masking of your desired object to the India domain. Suppose, you want to mask object 'x', and IMD is your gridded data available from IMD.

load('IMD.rda') x_masked<-regrid(x,is=IMD) map(x_masked)

I hope that will work. If not, just knock again with the error. Yes, one more thing, you may face the memory issue, so I would like to suggest you to use the IMD-1 Deg data.

Cheers!

[monamidutta]

Hi Suman Thank you for your reply. Actually I have already tried this but not working. max(lon(x)) and min(lon(x)) is less than max(lon(imd)) and min(lon(imd)) as a result I think interpolation is not working. Suppose x is 0.25 X 0.25 degree when I am applying regrid for 1 degree it is doing well however, it is not able to mask the undesired region with regrid command.

[brasmus]

Hi, I have masked Mozambique the way you asked. First, I created a polygon that had the shape of the country:

# This code was used to generate the Moz-object that produces the shape of Mozambique

Moz <- function() {
  data("geoborders")
  i1 <- (geoborders$x > 30) & (geoborders$x < 41) &
    (geoborders$y > -27) & (geoborders$y < -10.5)
  coast <- data.frame(x=geoborders$x[i1],y=geoborders$y[i1])
  borders <- attr(geoborders,'borders')
  i2 <- (borders$x > 30) & (borders$x < 41) & (borders$y > -27) & (borders$y < -10)
  i2[is.na(i2)] <- TRUE
  borders <- borders[i2,]

  i3 <- (borders$x < 31.7) & (borders$y < -25.6) |
    (borders$x < 31.3) & (borders$y > -25) & (borders$y < -20) |
    (borders$x < 30.35) & (borders$y > -16) & (borders$y < - 15.5) |
    (borders$x < 33.75) & (borders$y > -14) |
    (borders$x < 35) & (borders$y > -11.5)
  i3[is.na(i3)] <- TRUE
  borders[i3,] <- NA
  moz <- rbind(as.data.frame(coast),as.data.frame(borders))
  ok <- is.finite(moz$x) & is.finite(moz$y)
  moz <- moz[ok,]

  ## The ordering of borders and coasts is a bit random. Try to get a single
  ## thread.
  Moz <- moz
  print(moz[1,])
  Moz[-1,] <- NA; moz[1,] <- NA
  n <- length(Moz[,1])
  D <- rep(NA,n); D[1] <- 0
  for (i in 2:n) {
    d <- (Moz$x[i -1] - moz$x)^2 + (Moz$y[i-1] - moz$y)^2
    if (sum(is.finite(d))>0) {
      ii <- (1:n)[(d == min(d,na.rm=TRUE))]
      ii <- ii[is.finite(ii)][1]
      if (length(ii)>0) {
        D[i] <- min(d,na.rm=TRUE)
        Moz[i,] <- moz[ii,]
        moz[ii,] <- NA
        #print(Moz[i,])
      }
    }
  }
  Moz <- Moz[D < 0.1,]
  n <- length(Moz[,1])
  Moz <- Moz[1:(n-2),]
  Moz <- rbind(Moz,Moz[1,])
  save(Moz,file='Moz.rda')
}

## Run Moz
Moz()

Once the shape was made, I could use it to make a map with the shape of Mozambique:

## Plot map of gridded precipitation in Mozambique
## Rasmus.Benestad  2019-02-05

library(esd)
library(readINAMdata)
data("Moz")
if (file.exists('ERA5_annual-tp.rda')) load('ERA5_annual-tp.rda') else
   tp <- annual(retrieve('ERA5_tp_Mozambique_1979-2018.nc'),FU
N='sum')
lons <- rep(lon(tp),length(lat(tp)))
lats <- sort(rep(lat(tp),length(lon(tp))))
osmoz <- !is.inside(data.frame(x=lons,y=lats),Moz)
plot(lons,lats)
points(lons[osmoz],lats[osmoz],pch=19,col='red')
points(lons[!osmoz],lats[!osmoz],pch=19,col='green')
points(lons[is.na(osmoz)],lats[is.na(osmoz)],pch=19,col='black')
lines(Moz)

tp[,osmoz] <- NA
tp <- tp*1000
attr(tp,'unit') <-'mm'

load('CHIRPS_annual_1981-2018.rda')
lons <- rep(lon(TP),length(lat(TP)))
lats <- sort(rep(lat(TP),length(lon(TP))))
osmoz <- !is.inside(data.frame(x=lons,y=lats),Moz)

TP[,osmoz] <- NA
TP <- TP
attr(TP,'unit') <-'mm'

## plot the maps
mtp <- map(subset(tp,it=c(1981,2018)),FUN='mean',
           colbar=list(breaks=seq(0,2400,by=200)))
text(35,-10,'ERA5 mean annual rainfall')
mTP <- map(annual(TP,FUN='sum'),FUN='mean',colbar=list(breaks=seq(0,2400,by=200)
))
text(35,-10,'CHIRPS mean annual rainfall')
dtp <- regrid(mtp,is=mTP) - mTP
map(dtp,colbar=list(pal='t2m',breaks=seq(-1000,1000,by=200),rev=TRUE))
text(35,-10,'ERA5 -CHIRPS rainfall difference')

You need to do some manual work by first fixing the borders so that the shape looks OK.