[R-sig-Geo] Problems with image2Grid

[Roger Bivand]

On Thu, 24 Sep 2009, Michael Sumner wrote:

> Hello, Roger's reply to this made me realize I replied only to
> fernando, so I wanted to correct my first message and add some more.
>
> This reminds me that image2Grid needs check for irregular x or y
> values - it assumes the relationship of the first 2 in each define
> them all, so it should error in this case.

Sorry, I was looking at the x's - the longitude, not the y's where the 
problem was. Had the commands been included, it would have been easier to 
deconstruct. I'm committing a test in image2Grid() which gains a digits= 
argument to set the tolerance of uniqueness of step difference testing. 
With the check in place, the example data fails.

Paul Hiemstra put an interpolation on the R-wiki yesterday at:

http://wiki.r-project.org/rwiki/doku.php?id=tips:spatial-data:change_crs

perhaps you could add your example to it as well as posting?

Roger

>
> I'll try to put together an example that doe the resampling in R with
> overlay, or its underlying functions, to avoid external calls to GDAL.
>
> The y values in the example data are irregular and so this is not supported by
> Spatial[Grid/Pixels]DataFrame - although this is supported by image()
> for xyz lists of this type.
>
> range(diff(met$y))
> [1] 0.09527307 0.13665998
>
> You can convert this to points like this:
>
> library(sp)
> x <- data.frame(expand.grid(x = met$x, y = met$y), z = as.vector(met$z))
> coordinates(x) <- ~x+y
>
> If you know the projection and there is a regular transform that works
> - this is common in NetCDF files and environmental models where an
> underlying Mercator grid is use in lat/long form - only Y is
> "irregular". Having a wild guess:
>
> library(rgdal)
> proj4string(x) <- CRS("+proj=longlat +ellps=WGS84")
> x2 <- spTransform(x, CRS("+proj=merc +ellps=WGS84"))
> gridded(x2) <- TRUE
>
> image(x2)
>
> That works to create a regular grid and you could transform companion
> datasets in the same way, but may not be useful depending on your
> application. If you need to warp the grid to regular raster you can do
> that externally - with the GDAL command line utilities, which I have
> in my path. E.g.
>
> writeGDAL(x2, "file2convert.tif")
> system("gdalwarp file2convert.tif LL.tif -t_srs \"+proj=longlat +ellps=WGS84\"")
>
> ll.x <- readGDAL("LL.tif")
> summary(ll.x)
> Object of class SpatialGridDataFrame
> Coordinates:
>       min       max
> x -77.35186 -64.60699
> y -58.38184 -41.03831
> Is projected: FALSE
> proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
> Number of points: 2
> Grid attributes:
> cellcentre.offset  cellsize cells.dim
> x         -77.28616 0.1313904        97
> y         -58.31615 0.1313904       132
> Data attributes:
>   Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
> 0.003633 0.438200 0.438200 0.426500 0.438200 0.983700
>
> image(ll.x)
> contour(met, add = TRUE)
>
> HTH
>
> Regards, Mike.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no