[R-sig-Geo] how to calculate with grids?

[Robert J. Hijmans]

Bastian,

These rasters have a different extent and origin. So you cannot just
subtract them. This often some unfortunate processing choices earlier
on. But if that is not the case, I think you need to resample (and
hence distort) your data. You can use the raster package for this. Had
the origin been the same you could have used the 'crop' function only,
and not 'resample'.

install.packages("raster", repos="http://R-Forge.R-project.org")
library(raster)

dem1 <- raster("grid1.asc")
dem2 <- raster("dem2.asc")
intersect <- crop(dem2, dem1)
dem1r <- resample(dem1, intersect, method='bilinear', progress='text')
dem3 <- intersect - dem1r
plot(dem3)

#or vice versa, resample to dem2 rather than to dem1
intersect2 <- crop(dem1, dem2)
dem2r <- resample(dem2, intersect2, method='bilinear', progress='text')
dem4 <- dem2r - intersect2
plot(dem4)

Good luck,
Robert

On Tue, Mar 16, 2010 at 12:49 PM, Bastian Pöschl <rotate at gmx.li> wrote:
> Dear useRs,
>
>        I try to subtract one ascii grid from another.
>
>        By now i havn't managed to get a proper output.
>        Sorry I am not able to provide the files for testing because
>        they are very big.
>
>        The grids are read by
>
>        dem1<-readAsciiGrid("grid1.asc", as.image = FALSE, plot.image =
>        FALSE, colname = "height_old", proj4string = crs)
>        dem2<-readAsciiGrid("dem2.asc", as.image = FALSE, plot.image =
>        FALSE, colname = "height_new", proj4string = crs)
>
>        The cell dimension is not equal (i tried to convert them into
>        SpatialPointDataFrames to subtract the data, but struggled
>        because of the different extensions.
>        The cellsize is equal
>
>        1> summary(dem1)
>        Object of class SpatialGridDataFrame
>        Coordinates:
>                      min     max
>        coords.x1 3493070 3495230
>        coords.x2 5404550 5405920
>        Is projected: TRUE
>        proj4string :
>        [+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
>        +ellps=bessel +datum=potsdam +units=m +no_defs]
>        Number of points: 2
>        Grid attributes:
>          cellcentre.offset cellsize cells.dim
>        1         3493070.5        1      2160
>        2         5404550.5        1      1370
>        Data attributes:
>              Min.    1st Qu.     Median       Mean    3rd Qu.
>        Max.       NA's
>          420.0000   433.4856   443.1912   442.5298   450.7338
>        465.0000 80438.0000
>
>        1> summary(dem2)
>        Object of class SpatialGridDataFrame
>        Coordinates:
>                        min       max
>        coords.x1 3492999.5 3495000.5
>        coords.x2 5403999.5 5406000.5
>        Is projected: TRUE
>        proj4string :
>        [+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
>        +ellps=bessel +datum=potsdam +units=m +no_defs]
>        Number of points: 2
>        Grid attributes:
>          cellcentre.offset cellsize cells.dim
>        1           3493000        1      2001
>        2           5404000        1      2001
>        Data attributes:
>            Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
>        396.3400 425.7900 439.0900 437.6107 449.2700 470.1700
>
>
>
>        I found some hints in the help archive.
>        http://finzi.psych.upenn.edu/Rhelp10/2009-February/187216.html
>        tells about combining grids in an array, what i tried, but
>        failed.
>
>        I hope somebody can help me and give some advice how grids can
>        be merged and/or calculations can be don within the grid format.
>
>        This would be great.
>
>        Thank you,
>
>        Bastian
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>