I'm trying to use R to rasterize a number of shapefiles and I'm having trouble dealing with cells on the edge of shapes, in other words the if a cell is partially covered by an NA value and partially covered by a value associated with the shape. If any part of the cell is covered by a shape I would like the final raster product to have that value. For example in the image below I show the original shapefile (in green colors) and the result of a rasterize operation using R (in grey-scale). All of the red circles show areas where I would like the raster to have a value because some portion of the cell is covered by a shape.

enter image description here

I've tried using various functions with rasterize (max, mean, na.rm = T) and haven't had any luck producing the result I'm looking for. Below are links to the datafiles and my basic procedure.

raster: ftp://ftp.cdc.noaa.gov/Datasets/cpc_us_precip/precip.V1.0.mon.ltm.nc shapeFile: http://droughtmonitor.unl.edu/data/shapefiles_m/USDM_20110802_M.zip

R is my tool of choice but I'm happy with solutions in Arc or QGis as well. With the caveat that any solution needs to be relatively fast (I have 900 of these and I would ideally like to get this done in less than a day) and batchable (I'm familiar with model builder in Arc)

##Load libraries
library("rgdal") #Handling Shapefiles
library("raster") #Handling Raster files
library("ncdf4") #For reading raster from .nc file

##Load in shapefile as spatialpolygon
myShape <- readOGR(dsn = "USDM_20110802.shp", layer = "USDM_20110802")

##Load in Raster
#This will throw a waring about the variable name it can be ignored
#The netCDF file is multilayer but we only need on for a mask
myRaster <- raster("precip.V1.0.mon.ltm.nc")

##Rotate Raster-orignial file use longitude of 0 to 360 we want -180 to 180 to work with the shapefile
#This will also throw an error which can be ignored
myRaster <- rotate(myRaster)

##Rasterize the shape file
#I've tried all the following
shapeRaster <- rasterize(myShape, myRaster, field = "DM")
shapeRaster <- rasterize(myShape, myRaster, field = "DM", na.rm = T)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = max)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = max, na.rm = T)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = mean)

##My goal is for this cell to have a value of 0 instead of NA 
#in the image it's the top circled cell on the lefthand side
  • This is from the documentation from R and maybe a work around "For polygons, values are transferred if the polygon covers the center of a raster cell. For lines, values are transferred to all cells that are touched by a line. You can combine this behaviour by rasterizing polygons as lines first and then as polygons." Have you tried the Polygon to Raster function/tool in Arc? I'm not sure how it handles these cells and I can't find anything in the documentation. But it would be "easy" to process all 900 of them in Arc using Python/Arcpy. – F_Kellner Apr 27 '16 at 16:26
  • I actually had the same thought right after I posted this. It kinda worked with a few problems, the first is its slow...for some reason R rasterizes a polygon in a few seconds but took about 5 minutes to do the polyline. Second it creates weird behavior where polygons touch each other. Most of these files have polygons within polygons and the line to raster essentially expands all of the inner polygons by one grid cell. I think that makes sense. I've tried Arc and it seems limited to maximum area and cell center. – Adam Mccurdy Apr 27 '16 at 16:49
  • Would increasing the cell size grab up some of those edge polygons? – F_Kellner Apr 27 '16 at 17:37

Your Answer

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

Browse other questions tagged or ask your own question.