[R-sig-Geo] raster / polygon overlay

Fri Jan 8 14:58:33 CET 2010

Dear Laure,

If this is still actual, here are few tips on how to aggregate continuous information (rasters) given some polygon maps.

---------------------
# obtain a polygon map e.g. World countries:
> load(url("http://spatial-analyst.net/DATA/WorldPolyCountries.Rdata"))
> str(WorldPolyCountries at data)
'data.frame':   246 obs. of  11 variables:
 $ FIPS     : Factor w/ 243 levels "AA","AC","AE",..: 2 5 6 7 8 10 11 12 13 17 ...
 $ ISO2     : Factor w/ 246 levels "AD","AE","AF",..: 4 61 17 6 7 9 12 11 14 24 ...
 $ ISO3     : Factor w/ 246 levels "ABW","AFG","AGO",..: 15 64 18 6 11 3 12 10 16 25 ...
 $ UN       : int  28 12 31 8 51 24 16 32 36 48 ...
 $ NAME     : Factor w/ 246 levels "Afghanistan",..: 10 4 16 3 12 7 5 11 14 18 ...
 $ AREA     : int  44 238174 8260 2740 2820 124670 20 273669 768230 71 ...
 $ POP2005  : int  83039 32854159 8352021 3153731 3017661 16095214 64051 38747148 20310208 724788 ...
 $ REGION   : int  19 2 142 150 142 2 9 19 9 142 ...
 $ SUBREGION: int  29 15 145 39 145 17 61 5 53 145 ...
 $ LON      : num  -61.78 2.63 47.4 20.07 44.56 ...
 $ LAT      : num  17.1 28.2 40.4 41.1 40.5 ...
> writeOGR(WorldPolyCountries, "WorldPolyCountries.shp", "WorldPolyCountries", "ESRI Shapefile")

# download some rasters, e.g. global elevations:
> download.file("http://spatial-analyst.net/worldmaps/globedem.zip", destfile=paste(getwd(), "/glodedem.zip", sep=""))
> unzip(paste(getwd(), "/glodedem.zip", sep=""))
> rsaga.esri.to.sgrd(in.grids="globedem.asc", out.sgrd="globedem.sgrd", in.path=getwd())

# Get grid Statistics for Polygons using SAGA GIS:
# rsaga.get.usage("shapes_grid", 2)
> rsaga.geoprocessor(lib="shapes_grid", module=2, param=list(GRID="globedem.sgrd", POLY="WorldPolyCountries.shp", RESULT="WorldPolyCountries.shp", QUANTILES=TRUE, QUANTILE_STEP=2))

# To have more control, read maps R:
> worldgrid <- readGDAL("globedem.asc")
> names(worldgrid) <- "globedem"
> pixelsize <- worldgrid at grid@cellsize[1]

# convert Polygons to rasters:
> rsaga.geoprocessor(lib="grid_gridding", module=0, param=list(GRID="WorldCountries.sgrd", INPUT="WorldPolyCountries.shp", FIELD=3, LINE_TYPE=0, TARGET_TYPE=0, USER_CELL_SIZE=pixelsize, USER_X_EXTENT_MIN=worldgrid at bbox[1,1]+pixelsize/2, USER_X_EXTENT_MAX=worldgrid at bbox[1,2], USER_Y_EXTENT_MIN=worldgrid at bbox[2,1]+pixelsize/2, USER_Y_EXTENT_MAX=worldgrid at bbox[2,2]-pixelsize/2))
# read to R:
> worldgrid$UN <- as.factor(readGDAL("WorldCountries.sdat")$band1)

# aggregate:
> globedem.stats <- boxplot(worldgrid$globedem ~ worldgrid$UN, plot=FALSE)
# match the country names:
> UN.name <- merge(data.frame(UN=as.integer(globedem.stats$names)), WorldPolyCountries at data[,c("UN","NAME")], by="UN")
> globedem.stats$NAME <- UN.name$NAME
> data.frame(elev=globedem.stats$stats[3,], country=globedem.stats$NAME)[order(globedem.stats$stats[3,], decreasing=TRUE)[1:10],]
    elev     country
204 3396  Tajikistan
19  2802      Bhutan
114 2738  Kyrgyzstan
117 2192     Lesotho
6   2082     Andorra
16  1876     Armenia
85  1774   Greenland
174 1620      Rwanda
1   1588 Afghanistan
142 1566       Nepal

# 10 highest countries in the World...
---------------------

If you need to downscale gridded maps, then consider also using SAGA GIS, which has a very flexible down/up-scaling functionality:

http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE#Rescaling_maps_in_SAGA_GIS 

HTH

T. Hengl
http://home.medewerker.uva.nl/t.hengl/ 

> -----Original Message-----
> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-
> bounces at stat.math.ethz.ch] On Behalf Of laure velez
> Sent: Tuesday, January 05, 2010 2:47 PM
> To: r-sig-geo at stat.math.ethz.ch
> Subject: Re: [R-sig-Geo] raster / polygon overlay
> 
> 
> 
> Thanks again for the fast reply,
> 
> 
> 
> In fact this problem is only a tiny part of the project and I'll need
> 
> to run this overlay operation over more than 600 polygons (watersheds)
> 
> to get mean values (and standard error) of ~10 environmental variables
> 
> (rasters that are not at the same resolution, i.e 2.5 degrees or 0.1
> 
> degrees). Furthermore I will have to overlay the watershed polygons
> 
> with polygons of landcover.
> 
> 
> 
> The scale of this project is regional (the whole Mediterranean basin)
> 
> thus the number (and size) of rasters preclude analyses to be run at
> 
> the entire scale of the project.
> 
> 
> 
> Thus I can imagine the following automated process :
> 
> 
> 
> 1. calculate a buffer around the watershed polygon.
> 
> 
> 
> 2. subset the grid area that overlay with the buffer (will it work
> 
> with big grid cells ?)
> 
> 
> 
> 3. spsample the raster to a finer resolution (the number of new
> 
> samples will depend on the initial resolution of the raster), what
> 
> about using the bb argument of the spsample function ?
> 
> 
> 
> 3. overlay this resampled raster with the watershed polygon.
> 
> 
> 
> 4. get statistics (mean, sd, ...) from the resulting dataset.
> 
> 
> 
> As a newbie I would like to get some help (function names) for those
> 
> different steps (specifically 1 (buffer), 2 (grid subset)).
> 
> 
> 
> Thanks again for all your useful advices,
> 
> 
> 
> Laure.
> 
> _________________________________________________________________
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo