5

I am working on mapping some netcdf data in R and have been having a tough time trying to perform a clip with an imported country shapefile that I'm using. Here is my question from yesterday. If you know anything about that, feel free to pop over there and answer that question. :)

BUT I'll start fresh here. 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. The data is in netcdf format. I was able to map the netcdf data with the levelplot function from the rasterVis package. I was also able to map the netcdf data with a simple image function. Here's my code (sorry for the funky naming conventions; I've just been playing around with this):

library(maptools)
library(ncdf4)
library(RColorBrewer)
library(raster)

setwd("D:\\stuff")

#Madagascar shapefile
madagascar <- readShapeSpatial("MDG_adm0.shp")
plot(madagascar)

#Call on netcdf file
ncdf.data <- nc_open("swio_rpmaps_200_83.nc")

#Call on the 10 year wind return periods
returns <- ncvar_get(ncdf.data, "y010")

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

#plot image of windspeeds
image(lon, lat, returns, col = cm.colors(9,alpha=.6), add = TRUE)
#plot Madagascar on top of windspeeds
plot(madagascar, add = T)

Here is the result image: enter image description here

Is there any way to clip the Madagascar shapefile to this image? Or another way to manipulate the data so that I can perform a clip? Ideally, I would only have the data associated within the bounds of whatever country I'm mapping.

  • 1
    I think you're looking for the raster::mask function. – hrbrmstr Oct 21 '15 at 1:33
  • Any more details on how I should go about doing this? I'm very new to R and the help menu can be less than helpful most times. I'm playing around with it with no luck. The image doesn't seem to be inherently spatial, and the shapefile is defined as a SpatialPolygonsDataFrame. – Lauren Oct 21 '15 at 16:02
3

Picking up on the suggestion by @hrbrmstr here is the code. Since you you don't provide a sample file, I try to mimic your situation and create a sample the way you might be able to do it with your netcdf file. The shapefile I extracted from Natural Earth.

library(maptools)
library(raster)
madagascar <- readShapeSpatial("MDG_adm0.shp")

## make up some values 
# I am trying to replicate the way you would do it with the netcdf file..
lon <- seq(42, 51, .05) 
lat <- seq(-29, -10, .1)
returns <- runif(length(lat) * length(lon))

image(lon, lat, matrix(returns, length(lon), length(lat)), col = cm.colors(9,alpha=.6))
plot(madagascar, add = T)

enter image description here

Next we use the values from the netcdf file to create an actual raster object (a RasterLayer to be exact) that you can work with. To do this we need to set the dimensions of the raster (nrows and ncols) as well as the min and max values for both dimensions (x and y). We then use the setValues function to add the values (in your case windspeeds) for each raster cell.

UPDATE: Instead of doing this, create a Raster object directly as outlined by @RobertH .

# create the raster object using the values from above
r <- raster(nrows=length(lon), ncols=length(lat), 
        xmn=min(lon), xmx=max(lon),
        ymn=min(lat) ,ymx=max(lat))
r <- setValues(r, returns)

plot(r)
plot(madagascar, add = T)

enter image description here

Lastly, we use the mask function to clip the raster. This function takes any Spatial* object as mask, a SpatialPolygonsDataframe (like your madagascar one) being one of them.

#clip
r.clipped <- mask(r, madagascar)
plot(r.clipped)

enter image description here

  • thank you thank you thank you! This works. Would you mind explaining a bit what you did here? I'm an R newbie. However, it flips all of my "background" raster data 90 degrees to the right. Any idea what's going on there? I'd post a picture, but unsure how to in this comment window. – Lauren Oct 21 '15 at 19:06
  • @Lauren: I added more explanations into the post. Re your flipped image, make sure that you have lat and lon correctly assigned. Remember that x = lon and y = lat. You cannot post a picture in a comment, but you can post a link to a picture, if you upload it somewhere else. If you still have trouble, consider a new post and consider providing a sample of your netcdf file with your code. It makes it much easier for others to troubleshoot. – cengel Oct 21 '15 at 21:07
  • Thank you again! As far as the flipped image, I'm still perplexed. The first image call posts the correct orientation of the netcdf data. As soon as I rasterize it with your code, it flips 90 degrees. Very confusing! Would you recommend I post another question? – Lauren Oct 22 '15 at 15:03
  • Yes. And do provide a sample of your netcdf file. Perhaps the cell values are not assigned correctly. Also try what RobertH suggests. – cengel Oct 22 '15 at 22:33
  • There is still a lot of white space on both left and right of the cropped image. How to reduce/remove it? – Asad Ali Feb 1 '18 at 10:30
1

Adding to cengel's answer, to show how to get access to ncdf file as a RasterBrick object, and then contine:

library(raster)
madagascar <- shapefile("MDG_adm0.shp")
b <- raster("swio_rpmaps_200_83.nc", var= "y010")

bb <- crop(b, madagascar)
bb <- mask(bb, madagascar)
  • That's even better, as it bypasses my steps of manual extraction from the netcdf. Thank you for pointing this out! – cengel Oct 22 '15 at 22:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.