And as a tasty GRASS-based follow-up:
the interested user can take a look at the following functions;
# change the resolution, resampling is automatic via NN
g.region
r.mapcalc "new_resampled = old_map"
# up-sample (or down-scale) via NN, bilinear, or cubic interpolation
r.resamp.interp
# down-sample (or up-scale) via block- mean, median, mode, etc.
r.resamp.stats
# coordinate transformation across projections / datums
# will also resample to new grid topology via NN, bilinear, cubic,
# cubic-splines, etc.
r.proj
These operations can be performed from within R, via the spgrass6 package.
FOSS GIS means lots of options for these kind of common operations on
geographic data!
Cheers,
Dylan
On Monday 27 October 2008, Hengl, T. wrote:
> 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 at bbox[2,2]+meuse.grid.ll at bbox[2,1])/2*pi/180))/2 #
> > correction factor based on the latitude geogrd.cell <-
> > corrf*(meuse.grid.ll at bbox[1,2]-meuse.grid.ll at bbox[1,1])/meuse.grid at grid@c
> >ells.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 at coords[geoarc at coords[,2]==geoarc at coords[1,2],1]
> > grid.y <- geoarc at coords[geoarc at coords[,1]==geoarc at coords[1,1],2]
>
> # bilinear resampling:
> > library(akima)
> > dist.grid.ll <- interp(x=meuse.grid.ll at coords[,1],
> > y=meuse.grid.ll at 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 at 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 at uni-muenster.de]
> Sent: Fri 24/10/2008 5:24 PM
> To: Hengl, T.
> Cc: 'Pieter Beck'; r-sig-geo at 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 at grid@cells.dim[[2]],"
> >
> > ",geoarc at grid@cells.dim[[1]]," -crdsys=LatlonWGS84
> > -lowleft=(",geoarc at grid@cellcentre.offset[[1]],",",geoarc at grid@cellcentre
> >.offset[[2]],") -pixsize=",geoarc at 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 at stat.math.ethz.ch
> > [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Pieter Beck
> > Sent: Thursday, October 23, 2008 9:15 PM
> > To: r-sig-geo at 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 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
--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341