[GIS] Is it possible to clip a shapefile to an image() in R

clippolygonrshapefile

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:

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.

Best Answer

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)

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)

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)

Related Solutions

R Performance Polygons – Speeding Up Plotting of Polygons in R

I found 3 ways to increase the speed of plotting the country borders from shape files for R. I found some inspiration and code from here and here.

(1) We can extract the coordinates from the shape files to get the longitude and latitudes of the polygons. Then we can put them into a data frame with the first column containing the longitudes and the second column containing latitudes. The different shapes are separated by NAs.

(2) We can remove some polygons from our shape file. The shape file is very, very detailed, but some of the shapes are tiny islands that are unimportant (for my plots, anyways). We can set a minimum polygon area threshold to keep the bigger polygons.

(3) We can simplify the geometry of our shapes using the Douglas-Peuker algorithm. The edges of our polygon shapes can be simplified, as they are very intricate in the original file. Fortunately, there is a package, rgeos, that implements this.

Set up:

# Load packages
library(rgdal)
library(raster)
library(sp)
library(rgeos)

# Load the shape files
can<-getData('GADM', country="CAN", level=0)
usa<-getData('GADM', country="USA", level=0)
mex<-getData('GADM', country="MEX", level=0)

Method 1: Extract the coordinates from the the shape files into a data frame and plot lines

The major disadvantage is that we lose some information here when compared to keeping the object as a SpatialPolygonsDataFrame object, such as the projection. However, we can turn it back into an sp object and add back the projection information, and it is still faster than plotting the original data.

Note that this code runs very slowly on the original file because there are a lot of shapes, and the resulting data frame is ~2 million rows long.

Code:

# Convert the polygons into data frames so we can make lines
poly2df <- function(poly) {
  # Convert the polygons into data frames so we can make lines
  # Number of regions
  n_regions <- length(poly@polygons)
  
  # Get the coords into a data frame
  poly_df <- c()
  for(i in 1:n_regions) {
    # Number of polygons for first region
    n_poly <- length(poly@polygons[[i]]@Polygons)
    print(paste("There are",n_poly,"polygons"))
    # Create progress bar
    pb <- txtProgressBar(min = 0, max = n_poly, style = 3)
    for(j in 1:n_poly) {
      poly_df <- rbind(poly_df, NA, 
                       poly@polygons[[i]]@Polygons[[j]]@coords)
      # Update progress bar
      setTxtProgressBar(pb, j)
    }
    close(pb)
    print(paste("Finished region",i,"of",n_regions))
  }
  poly_df <- data.frame(poly_df)
  names(poly_df) <- c('lon','lat')
  return(poly_df)
}

Method 2: Remove small polygons

There are many small islands that are not very important. If you check some of the quantiles of the areas for the polygons, we see that many of them are miniscule. For the Canada plot, I went down from plotting over a thousand polygons to just hundreds of polygons.

Quantiles for the size of polygons for Canada:

          0%          25%          50%          75%         100% 
4.335000e-10 8.780845e-06 2.666822e-05 1.800103e-04 2.104909e+02

Code:

# Get the main polygons, will determine by area.
removeSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons, 
                  function(x) sapply(x@Polygons, function(y) y@area))
  
  # Quick summary of the areas
  print(quantile(unlist(areas)))
  
  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))
  
  # Get only the big polygons and extract them
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

Method 3: Simplify the geometry of the polygon shapes

We can reduce the number of vertices in our polygon shapes using the gSimplify function from the rgeos package

Code:

can <- getData('GADM', country="CAN", level=0)
can <- gSimplify(can, tol=0.01, topologyPreserve=TRUE)

Some benchmarks:

I used elapsed system.time to benchmark my plotting times. Note that these are just the times for plotting the countries, without the contour lines and other extra things. For the sp objects, I just used the plot function. For the data frame objects, I used the plot function with type='l' and the lines function.

Plotting the original Canada, USA, Mexico polygons:

73.009 seconds

Using Method 1:

2.449 seconds

Using Method 2:

17.660 seconds

Using Method 3:

16.695 seconds

Using Method 2 + 1:

1.729 seconds

Using Method 2 + 3:

0.445 seconds

Using Method 2 + 3 + 1:

0.172 seconds

Other remarks:

It seems that the combination of methods 2 + 3 gives sufficient speed ups to the plotting of polygons. Using methods 2 + 3 + 1 adds the problem of losing the nice properties of sp objects, and my main difficulty is applying projections. I hacked something together to project a data frame object, but it runs rather slow. I think using method 2 + 3 provides sufficient speed ups for me until I can get the kinks out of using method 2 + 3 + 1.

[GIS] R – Plotting shapefile with geographical (raster) data

To me there are really 2 unrelated questions here.

A) What is the projection of my shapefile and how can I merge/plot it with another layer.

B) What is the form of these raster data, and why are there in a table when they could actually be a grid?

There's a lot of possibilities here and not enough to go on. I'd like to know the answers to these:

I would find out the projection of the shapefile, and you cannot do that with maptools::readShapeLines.

1) read it with rgdal and transform:

library(rgdal)
ned.lines <- readOGR("d:/readShapefiles", "NLD_adm/NLD_adm1")

## this will tell you what the projection is
## find out and tell us
 proj4string(ned.lines)

You might find this makes it work with your grid data, but there's a lot of guessing here:

plot(spTransform(ned.lines, CRS("+proj=longlat +datum=WGS84")))

2) Find out the actual projection and get back to us, you might need this if you cannot install rgdal.

readLines("d:/readShapefiles/NLD_adm/NLD_adm1.prj")

(Note that this is reading text lines, not geometry).

If that fails because the file doesn't exist, then there's only more guesses we can do. It looks like you don't have a grid of the radar data, so I would go back to the source and get an actual grid, or you'll need to rasterize the data from the irregular points.

Best Answer

Related Solutions

R Performance Polygons – Speeding Up Plotting of Polygons in R

[GIS] R – Plotting shapefile with geographical (raster) data

Related Question