¿Cómo recortar un polígono espacial a una gran matriz (que se traza como un gran enrejado, originado como un netcdf) en R?

Question

En primer lugar, me gustaría aclarar que soy un novato en SIG y que acabo de empezar a aprender R. Por lo tanto, si utilizo una terminología incorrecta o mi pregunta es confusa, le ruego que me ayude a aclararlo.

En mi trabajo tenemos previsto cartografiar toneladas y toneladas de datos meteorológicos de todo el mundo. Para practicar, la persona que nos suministra los datos nos ha proporcionado algunos datos "de ejemplo" con los que podemos jugar. A mi jefe le gustaría que yo creara mis mapas en R. Se pretende que estos mapas sean mapas de periodos de retorno ( https://en.wikipedia.org/wiki/Return_period ). Se supone que los datos están en formato netcdf.

He conseguido encontrar varios tutoriales que me han ayudado a mapear realmente los datos de netcdf. Aquí está el código:

#set working directory
setwd("D:\\stuff")

#add libraries
library(maptools)
library(ncdf4)
library(RColorBrewer)
library(lattice)
library(rasterVis)

#add madagascar shapefile
proj <- CRS('+proj=longlat +ellps=WGS84')
madagascar <- readShapeSpatial("MDG_adm0.shp", proj4string = proj)

#read netcdf data for 10 year return period
ncdf.data <- nc_open("swio_rpmaps_200_83.nc")
returns <- ncvar_get(ncdf.data, "y010")

#define latitude and longitude
lon <- ncdf.data$dim$longitude$vals
lat <- ncdf.data$dim$latitude$vals

grid <- expand.grid(lon = lon, lat = lat)
cutpts <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)
test <- levelplot(returns ~ lon * lat, 
          data = grid, 
          at = cutpts, 
          cuts = 10, 
          pretty = T, 
          col.regions = rev(brewer.pal(10, "RdBu")), 
          ylab = "Latitude",
          xlab = "Longitude",
          main = "10 Year Return Period")

#return periods with madagascar plotted on top
test + layer(sp.polygons(madagascar))

Este es el aspecto del mapa:

Problema: Quiero poder recortar SÓLO la parte de Madagascar de los datos de velocidad del viento presentados en este mapa. No puedo por mi vida averiguar cómo hacer esto. He intentado usar gIntersection del paquete rgeos, pero los datos netcdf (que aparecen como una gran matriz en mi espacio de trabajo; cuando los defino como "test" y uso levelplot, se convierten en un gran enrejado en el espacio de trabajo) aparentemente no son espaciales, así que no me deja hacer nada.

¿Hay alguna forma de realizar este clip con el código que he escrito? Si no, ¿cómo podría superponer un polígono espacial con estos datos netcdf? Esto me supera.

Answer 1

He aquí un planteamiento sencillo:

library(raster)
madagascar <- getData('GADM', country=MDG', level=0)
b <- brick("swio_rpmaps_200_83.nc", var= "y010")

bb <- crop(b, madagascar)
bb <- mask(bb, madagascar)

Para el trazado, véase

plot(bb, 1)
library(rasterVis)    
levelplot(bb)

He descargado el archivo y puede hacer que esto funcione mediante el uso de raster en lugar de brick (que falla con un archivo ncdf de una sola capa...)

r <- raster(f, varname='y010')
r
#class       : RasterLayer 
#dimensions  : 499, 699, 348801  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : 25.025, 59.975, -29.975, -5.025  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : swio_rpmaps_200_83.nc 
#names       : y010 
#zvar        : y010 

mdg <- getData('GADM', country='MDG', level=0)
x <- crop(r, mdg)
y <- mask(x, mdg)
plot(y)

Como has dicho que vas a hacer esto con frecuencia, yo crearía una máscara permanente para acelerar las cosas, y utilizaría una medida redondeada para evitar posibles problemas de alineación.

e <- floor(extent(mdg))
x <- crop(r, e)
x <- rasterize(mdg, x, filename='mdgmask.tif', datatype='INT2S')

Después de hacer eso una vez, ahora puede hacer

r <- brick(f, varname='y010')
m <- raster('mdgmask.tif')
x <- crop(r, m)
y <- mask(x, m)