Skip to content

calculate raster values based on vector regions

8 messages · Alexander.Herr at csiro.au, Tomislav Hengl, Kamran Safi +1 more

Original
Alexander.Herr at csiro.au
# 
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
Tomislav Hengl
# 
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:
# load the gridded map:
# load the shapefile:
# convert the polygon map to a raster map:
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))
# summary statistics per polygon class:
col=rainbow(length(levels(negotingrid$SOIL))))
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
Alexander.Herr at csiro.au
# 
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
Tomislav Hengl
# 
My appologies. The previous script had few (5) typos:
# load the gridded map:
# convert the polygon map to a raster map:
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))
# summary statistics per polygon class:
$polygons))))
HTH

Tom Hengl
http://spatial-analyst.net
Edzer Pebesma
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Alexander.Herr at csiro.au
    

    

      
      #
    

  


  

      
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: