[R-sig-Geo] GridTopology, Gluing SpatialGridDataFrames together from adjacent Rasters

[Robert J. Hijmans]

Chris,

Something along these lines might work:

library(raster)
drg85<- raster("~path/drg85.tif")
drg97<- raster("~path/drg97.tif")
drg85_p<- crop(drg85, extent(xmin, xmax, ymin, ymax))
drg97_p<- crop(drg97, extent(xmin2, xmax2, ymin2, ymax2))

ext <- unionExtent(drg85_p, drg97_p)
r1 <- expand(drg85_p, ext)
x <- crop(r1, drg97_p)
x <- resample(drg97_p, x, method='bilinear')

r <- merge(r1, y, filename='test.tif')


Robert




On Thu, Jul 22, 2010 at 10:03 AM, Chris English <sglish at hotmail.com> wrote:
>
> Hi List,
>
> I am new to R spatial.  I have group of islands in a lake that I kayak around and observe birds, animals, and vegetation.  I wanted to create a map of my trips and points of observation.
>
> My study area is essentially at the edges of two 7.5 quads, south edge of top and north edge of bottom.  My goal has been to join parts of these two quad maps and derive a map showing my study area.
>
> To Start:
>
> drg85<- raster("~path/drg85.tif")
> drg87<- raster("~path/drg97.tif")
>
> drg85_p<- crop(drg85, extent(xmin, xmax, ymin, ymax))
> if(require(rgdal)){drg85_pw<- writeRaster(drg85_p, filename="drg85_p",
> format="GTiff", overwrite=TRUE)}
>
> (do the same crop and write for drg97_p)
> the cropping also allowed for plotting in a small RAM laptop
> drg85_p<- readGDAL("~path/drg85_p.tif")#returns SpatialGridDataFrame
> drg97_p<- readGDAL("~path/drg97_p.tif")#returns SpatialGridDataFrame
>

a more direct approach

drg85<- raster("~path/drg85.tif")
drg97<- raster("~path/drg97.tif")

drg85_p<- crop(drg85, extent(xmin, xmax, ymin, ymax) , filename =
"drg85_p.tif") # filename can be omitted
drg97_p<- crop(drg97, extent(xmin, xmax, ymin, ymax) , filename = "drg85_p.tif")

drg85_p<- as(drg85, 'SpatialGridDataFrame')
drg97_p<- as(drg97, 'SpatialGridDataFrame')



> So I tried to cbind:
>  drg8597_ps<- cbind(drg85_p, drg97_p)
> Error in stop.ifnot.equal(x, grds[[1]]) : topology is not equal
>
>> getGridTopology(drg85_p)
>                             x                    y
> cellcentre.offset 1.011904e+06   4.556105e+05
> cellsize              8.252375e+00   8.252375e+00
> cells.dim            2.001000e+03   8.960000e+02
>> getGridTopology(drg97_p)
>                             x                  y
> cellcentre.offset 1.011901e+06   4.54005e+05
> cellsize              8.251850e+00  8.25185e+00
> cells.dim            2.013000e+03   2.00000e+02
>
> And indeed they're different, most especially between y cells.dim.
> Checked to see if the grids held data:
>> fullgrid(drg85_p)
> [1] TRUE
>> fullgrid(drg97_p)
> [1] TRUE
> So I thought, well most of the area is in drg85_p, why not make the topology of drg97_p the same as drg85_p and tried:
>
>> drg97_p<- GridTopology(cellcentre.offset(x=1.011904e+06,y=4.556105e+05),cellsize(x=8.252375e+00, y=8.253275e+00), cells.dim(x=2.001000e+03 y=8.960000e+02))
> Error in initialize(value, ...) :
>  could not find function "cellcentre.offset"
>
> OK- get the form right:
>> drg97_p<- GridTopology(c(1.011904e+06, 4.556105e+05),c(8.252375e+00, 9.252375e+00),c(2.0010000e+03,8.960000e+02))
>> fullgrid(drg97_p)
> [1] FALSE
>
> Oops, the data.  Well, I've made a grid. Go readGDAL to drg97_pd again.
>
> And now I can't figure out how to read out the data and write to my equal topology.
> Any help would be appreciated.
>
> Chris
>
>
>
>
>
>
>
> _________________________________________________________________
> The New Busy is not the old busy. Search, chat and e-mail from your inbox.
>
> N:WL:en-US:WM_HMP:042010_3
>        [[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
>