[R-sig-Geo] calculate raster values based on vector regions
Alexander.Herr at csiro.au
Alexander.Herr at csiro.au
Thu Jan 29 00:18:08 CET 2009
Thanks Edzer,
That would do it. Unfortunately I can't read the (ESRI) grid, it is too large:
raster1 has GDAL driver AIG
and has 73772 rows and 80264 columns
Error in array(dim = as.integer(c(rev(output.dim), length(band)))) :
'dim' specifies too large an array
I 'll have to cut into smaller tiles before processing, unless I can find a way of doing polygon overlay in package{raster} on rows
Cheers
Herry
-----Original Message-----
From: Edzer Pebesma [mailto:edzer.pebesma at uni-muenster.de]
Sent: Thursday, January 29, 2009 2:38 AM
To: Herr, Alexander Herr - Herry (CSE, Gungahlin)
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] calculate raster values based on vector regions
Alexander, below is a reproducible example that uses a polygon data set
from package maptools. You should read your polygons e.g. through
readOGR (package rgdal) and your grid through readGDAL.
If your grid comes from readGDAL, then 2 and 3 can be omitted, and 5
simplifies to:
out$SID79 = nc$SID79[idx]
All in R. Hth,
--
Edzer
library(maptools)
# 1 read a SpatialPolygonsDataFrame
nc <- readShapePoly(system.file("shapes/sids.shp",
package="maptools")[1], proj4string=CRS("+proj=longlat
+datum=NAD27"))
# 2 sample points on a regular grid:
grd = spsample(nc, 5000, "regular", offset = c(.5,.5))
# 3 convert points into grid:
gridded(grd) = TRUE
# 4 find index of polygons on grid cell centres:
idx = overlay(grd, nc)
# 5 merge grid points with attributes from polygons:
out = SpatialPixelsDataFrame(grd, data.frame(SIDS79 = nc$SID79[idx]))
# 6 plot:
spplot(out, col.regions=bpy.colors())
Alexander.Herr at csiro.au wrote:
> Hi list,
>
> I try to explain my problem a bit better.
>
> 1) Vector and raster have the same extent.
>
> 2) Vector data:
> Several polygones with different attribute values
> 3) raster data:
> A raster with points over the areas (and also NAs for areas where it is not possible to have a value from polygones). These points need to get the polygon attribute values assigned
>
> Vector Raster after assiging polygon values to raster
> ----------- ------------ ------------
> | 3 | 1 | | +na + + na | |3na 3 1 na |
> |-----------| | na na na na| |na na na na|
> | 7 | 5 | |++na na + | |7 7na na 5 |
> ----------- ------------- -------------
>
> + denotes a raster value
>
> Does this make more sense?
>
> In gdal/fwtools is a function rgdal_rasterize, which could potentially do this, but I'd rather to all in R.
>
> Cheers and thanks
> Herry
>
> ________________________________________
> From: Tomislav Hengl [T.Hengl at uva.nl]
> Sent: Wednesday, January 28, 2009 7:56 PM
> To: r-sig-geo at stat.math.ethz.ch
> Cc: Herr, Alexander Herr - Herry (CSE, Gungahlin)
> Subject: RE: [R-sig-Geo] calculate raster values based on vector regions
>
> Dear Herry,
>
> If I understand what you problem, one solution is to use R+SAGA. You should first convert the
> polygon map to the same grid, and then you can load it to R and do any type of aggregation:
>
>
>> library(maptools)
>> library(rgdal)
>> library(RSAGA)
>>
> # load the gridded map:
>
>> rastermap <- readGDAL("rastermap.asc")
>>
> # load the shapefile:
>
>> rsaga.esri.to.sgrd(in.grids=" rastermap.asc", out.sgrds=" rastermap.sgrd", in.path=getwd())
>>
> # convert the polygon map to a raster map:
>
>> cellsize <- rastermap at grid@cellsize[1]
>> rsaga.geoprocessor(lib="grid_gridding", module=3, param=list(GRID="polygons.sgrd",
>>
> INPUT="polygons.shp", FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=cellsize,
> USER_X_EXTENT_MIN=rastermap at bbox[1,1]+cellsize, USER_X_EXTENT_MAX=rastermap at bbox[1,2]-cellsize,
> USER_Y_EXTENT_MIN=rastermap at bbox[2,1]+cellsize, USER_Y_EXTENT_MAX=rastermap at bbox[2,2]-cellsize))
>
>> rsaga.sgrd.to.esri(in.sgrds="polygons.sgrd", out.grids="polygons.asc", out.path=getwd(), prec=0)
>> rastermap$polygons <- as.factor(readGDAL("polygons.asc"))
>>
> # summary statistics per polygon class:
>
>> raster.polygons <- boxplot(band1 ~ polygons, rastermap @data,
>>
> col=rainbow(length(levels(negotingrid$SOIL))))
>
>> str(raster.polygons)
>>
>
> You will of course need to adjust the "FIELD" position in the attribute table and the precision
> "prec".
>
> To obtain and use SAGA, check this article:
> http://spatial-analyst.net/wiki/index.php?title=Software
>
> cheers,
>
> How is the weather now in Canberra? In Amsterdam is freezing brrr.
>
> 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 Alexander.Herr at csiro.au
>> Sent: Wednesday, January 28, 2009 2:10 AM
>> To: r-sig-geo at stat.math.ethz.ch
>> Subject: [R-sig-Geo] calculate raster values based on vector regions
>>
>>
>> Hi List,
>>
>> How would I go about assigning values to a raster based on polygon regions in R. Examples would be
>> most appreciated.
>>
>> I have (vector) regions assigned to a specific value. The raster has NAs and some pixels where
>> these values are likely to occur on ground. I need to assign these values to the raster and can do
>> this in ArcInfo through vectors and rasterizing but only to a limited raster size. And R is much
>> more preferable anyway...
>>
>> Any help appreciated.
>> Cheers
>> Herry
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
More information about the R-sig-Geo
mailing list