Problems with image2Grid

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.

--
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