Edzer,

I agree, your example is straight forward. But this is really a simplification of the problem, especially if you go from Cartesian to geographical coordinates (LatLon).

This is how typically reproject+resample grids in R:

> data(meuse.grid)
> coordinates(meuse.grid) <- ~x+y
> gridded(meuse.grid) <- TRUE
> proj4string(meuse.grid) = CRS("+init=epsg:28992")

# reprojection (first from grid to point map)
> meuse.grid.ll <- spTransform(meuse.grid, CRS("+proj=longlat +datum=WGS84"))

# estimate a suitable grid cell size in the new coordinate system:
> corrf <- (1+cos((meuse.grid.ll@bbox[2,2]+meuse.grid.ll@bbox[2,1])/2*pi/180))/2 # correction factor based on the latitude
> geogrd.cell <- corrf*(meuse.grid.ll@bbox[1,2]-meuse.grid.ll@bbox[1,1])/meuse.grid@grid@cells.dim[1]
> geoarc <- spsample(meuse.grid.ll, type="regular", cellsize=c(geogrd.cell,geogrd.cell))
> gridded(geoarc) = TRUE
> gridparameters(geoarc)
   cellcentre.offset     cellsize cells.dim
x1          5.721603 0.0004603493        96
x2         50.956995 0.0004603493        80
> gridparameters(meuse.grid)
> grid.x <- geoarc@coords[geoarc@coords[,2]==geoarc@coords[1,2],1]
> grid.y <- geoarc@coords[geoarc@coords[,1]==geoarc@coords[1,1],2]
# bilinear resampling:
> library(akima)
> dist.grid.ll <- interp(x=meuse.grid.ll@coords[,1], y=meuse.grid.ll@coords[,2], z=meuse.grid.ll$dist, xo=grid.x, yo=grid.y, extrap=F, linear=T)
> dist.grid.ll <- as.SpatialGridDataFrame.im(as.im(dist.grid.ll))
> proj4string(dist.grid.ll) = CRS("+proj=longlat +datum=WGS84")
> print(spplot(meuse.grid["dist"], scales=list(draw=T, cex=.8)), split=c(1,1,2,1), more=T)
> print(spplot(dist.grid.ll[1], scales=list(draw=T, cex=.8)), split=c(2,1,2,1), more=T)

What 'interp' method of akima does is exactly what e.g. MapResample method in ILWIS does (except ILWIS does it in one line, but you have to get the data to ILWIS format first).

I just remembered that one can also use SAGA GIS to reproject/resample grids from a coordinate system to geographic coordinates (and back) using the proj4 module:

> rsaga.get.modules("pj_proj4")
> rsaga.get.usage("pj_proj4",1)
...
SAGA CMD 2.0.3
library path:   C:/Progra~1/saga_vc/modules
library name:   pj_proj4
module name :   Proj4 (Grid)
author      :   (c) 2003 by O.Conrad

Note that SAGA makes it possible to use a range of resampling techniques (this would allow you to do e.g. down/up-scaling):

...
 -INTERPOLATION:<num>                  Grid Interpolation
        Choice
        Available Choices:
        [0] Nearest Neigbhor
        [1] Bilinear Interpolation
        [2] Inverse Distance Interpolation
        [3] Bicubic Spline Interpolation
        [4] B-Spline Interpolation

This is an example of a command that can be used to reproject/resample from TM to geographic grid (note that you can also leave it to SAGA to estimate the bounding box / new grid topology):

> rsaga.esri.to.sgrd(in.grids="dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
> rsaga.geoprocessor(lib="pj_proj4", module=1, param=list(SOURCE="dem25m_ll.sgrd", DIRECTION=1, PROJ_TYPE=103, ELLIPSOID=0, ELLPS_PREDEF=11, UNIT=1, LON_0=16.5, LAT_0=0, X_0=2500000, Y_0=0, OUT_GRID="dem25m_ll.sgrd", TARGET_TYPE=0, INTERPOLATION=1, GET_USER_SIZE=geoarc@grid@cellsize[[1]]))

# unfortunately this does not work because I am not able to define the correct geodetic datum (dX, dY, dZ) for baranja hill dataset (http://geomorphometry.org/data.asp) in SAGA, so I get an empty grid!

Unfortunately SAGA does not (yet!) support EPSG codes, so that the choice of coordinate systems is relatively limited.


Tom Hengl
http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE


-----Original Message-----
From: Edzer Pebesma [mailto:edzer.pebesma@uni-muenster.de]
Sent: Fri 24/10/2008 5:24 PM
To: Hengl, T.
Cc: 'Pieter Beck'; r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] subsetting
 
Tomislav Hengl wrote:
> In GIS terms, getting a raster from one grid topology to (any) other grid topology (also subgrids)
> is referred to 'resampling' (I do not think that this is yet implemented in any R package; a
> combination of akima and spTransform will do the trick but it is not really compact).
I do think so: overlay in sp is your friend:

library(sp)
data(meuse.grid)
gridded(meuse.grid) = ~x+y
grd = spsample(meuse.grid, n = 100, type = "regular")
image(meuse.grid["dist"])
points(grd, pch = 3)

gridded(grd)= TRUE
overlay(meuse.grid, grd)
grddf = SpatialPixelsDataFrame(grd, as(meuse.grid, 
"data.frame")[overlay(meuse.grid, grd),])
image(grddf["dist"])
points(as(grd, "SpatialPoints"), pch = 3)

Note that there is only a single code line that does the resample. If 
anyone has suggestions for a cleaner or simpler interface, or a function 
or method resample, please let me know.

This of course does not do any interpolation; I consider that a 
different chapter (and package). I cannot see how spTransform could help 
here.When warping grids, you'd have to revert them first into polygons 
and spTransform them, to conserve the non-square grid cells resulting.
--
Edzer


>  A generic
> function to combine R+ILWIS to put a map to a different grid (e.g. to LatLon coordinates) would be:
>
> # write to ILWIS:
>   
>> writeGDAL(dem25m[1], "dem25m.mpr", "ILWIS")
>>     
> # create a new grid:
> ...
>   
>> gridparameters(geoarc)
>>     
> # resample the map (Bilinear) to the new geographic grid:
>   
>> shell(cmd=paste(ILWIS, " crgrf geoarc.grf ",geoarc@grid@cells.dim[[2]],"
>>     
> ",geoarc@grid@cells.dim[[1]]," -crdsys=LatlonWGS84
> -lowleft=(",geoarc@grid@cellcentre.offset[[1]],",",geoarc@grid@cellcentre.offset[[2]],")
> -pixsize=",geoarc@grid@cellsize[[1]],sep=""), wait=F)
>   
>> shell(cmd=paste(ILWIS, "dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, BiLinear)"), wait=F)
>> shell(cmd=paste(ILWIS, "open dem25m_ll_c.mpr -noask"), wait=F)
>>     
> see http://geomorphometry.org/R.asp 
>
> hth,
>
> Tom Hengl
> http://spatial-analyst.net 
>
>
>
>
> -----Original Message-----
> From: r-sig-geo-bounces@stat.math.ethz.ch [mailto:r-sig-geo-bounces@stat.math.ethz.ch] On Behalf Of
> Pieter Beck
> Sent: Thursday, October 23, 2008 9:15 PM
> To: r-sig-geo@stat.math.ethz.ch
> Subject: [R-sig-Geo] subsetting
>
> Dear all,
>
> I am a new user to the world of spatial applications of R and am wondering
> about this trivial question:
>
> What is the easiest way to subset a SpatialGridDataFrame, based on 4 corner
> coordinates?
>
>  
>
> The reason I ask is that I want to do an overlay of a polygon file on a
> large raster dataset, to extract summary statistics from it. Both files are
> extremely large, however, and just applying overlay causes memory trouble. I
> was hoping to loop through the polygons, subset the SpatialGridDataFrame
> using the bounding box surrounding a polygon and then create the summary
> statistics one polygon at-a-time. To do this, I'd need to subset the
> SpatialGridDataFrame based on the bounding box corner coordinates.
>
> Thanks in advance for the help,
>
>  
>
> Kind Regards,
>
> Pieter Beck
>
>  
>
>
> 	[[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
>   

-- 
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@wwu.de




	[[alternative HTML version deleted]]


