[R-sig-Geo] Problems with image2Grid

[Michael Sumner]

I put it on the wiki page, hopefully not too obscure - I certainly
encounter that particular issue quite frequently.

Regards, Mike.

On Thu, Sep 24, 2009 at 8:37 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
>