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 <[email protected]> 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
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
