Friday, 7 October 2016

How to clip a spatial polygon to a large matrix (that is plotted as a large trellis, originated as a netcdf) in R?


I would first like to issue a disclaimer that I am a novice at GIS, and I have only just begun learning R. So if I use the wrong terminology for something or my question is confusing, please help clear that up for me.


So for my work, we are planning on mapping tons and tons of weather data all across the globe. For practice, the guy providing our data provided us with some "example" data to play with. My boss would like for me to create my maps in R. These maps are intended to be return period maps (https://en.wikipedia.org/wiki/Return_period). The data is supposed to be in netcdf format.


I've managed to find several tutorials that have helped me to actually map the netcdf data. Here's the code:


#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))

Here's what the map looks like:enter image description here


Problem: I want to be able to clip out JUST the Madagascar portion of the windspeed data presented in this map. I cannot for the life of me figure out how to do this. I attempted to use gIntersection from the rgeos package, but the netcdf data (which shows up as a large matrix in my workspace; when I define it as "test" and use levelplot, it becomes a large trellis in the workspace) is apparently not spatial, so it won't let me do anything.


Is there any way to perform this clip with the code I've written? If not, how would I be best suited to overlay a spatial polygon with this netcdf data? I'm in WAY over my head here!




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...