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