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(worldpolycountr...@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 <- worldg...@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=worldg...@bbox[1,1]+pixelsize/2,
> user_x_extent_max=worldg...@bbox[1,2],
> user_y_extent_min=worldg...@bbox[2,1]+pixelsize/2,
> user_y_extent_max=worldg...@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)),
> worldpolycountr...@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-boun...@stat.math.ethz.ch [mailto:r-sig-geo-
> boun...@stat.math.ethz.ch] On Behalf Of laure velez
> Sent: Tuesday, January 05, 2010 2:47 PM
> To: r-sig-geo@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@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
