0

Forgive me for me being naive. I am having trouble with re-projecting SGDF.

I have an xyz (x=longitude, y=latitude, z=value) 3-column dataset. Longitude and Latitude values are coordinates in EPSG:3035 format. I converted the data frame into a grid with a spatial resolution of 5 km * 5km projected in EPSG:3035. I wish to re-project the grid to EPSG:4326 with a spatial resolution of 0.05*0.05. However i get the following warning messages:

Warning messages:
1: In spTransform(radon, CRS("+init=epsg:4326")) :   Grid warping not available, coercing to points
2: In spTransform(as(x, "SpatialPixelsDataFrame"), CRSobj, ...) :   Grid warping not available, coercing to points

Can anybody please tell me, how i can re-project the grid. Below is a small reproducible example:

library(sp)
library(rgdal)
library(raster)

x=c(5013500, 5018500, 4883500, 4888500, 4893500, 4898500, 4908500,4948500, 4953500, 4958500, 4963500, 4973500, 4978500, 4988500, 5008500, 5013500, 5028500, 4878500, 4883500, 4888500, 4893500,4898500, 4903500,4928500, 4963500, 4968500, 4973500, 4978500, 4983500, 4988500)

y=c(5395500, 5395500, 5390500, 5390500, 5390500, 5390500, 5390500,5390500, 5390500, 5390500, 5390500, 5390500, 5390500, 5390500, 5390500,5390500, 5390500, 5385500, 5385500, 5385500, 5385500, 5385500, 5385500,5385500, 5385500, 5385500, 5385500, 5385500, 5385500, 5385500)

z=c(1.74, 1.74, 1.82, 1.82, 1.82, 1.81, 1.81, 1.78, 1.77, 1.77, 1.76,1.76, 1.75, 1.74, 1.73, 1.73, 1.72, 1.82, 1.82, 1.81, 1.81, 1.80, 1.80, 1.78, 1.75, 1.75, 1.74, 1.74, 1.73, 1.73)

df1=data.frame(x,y,z)
coordinates(df1) <- c("x", "y")
proj4string(df1)=CRS("+init=epsg:3035")
gridded(df1)=TRUE
fullgrid(df1)=TRUE

getGridTopology(df1)

                     x       y
cellcentre.offset 4878500 5385500
cellsize             5000    5000
cells.dim              31       3

newgrid = spTransform(df1, CRS("+init=epsg:4326"))

Well, thanks to Paul, i was able to re-project the grid using gdalwarp. But still, the spatial resolution is different:

    Coordinates:
    min      max
x -31.34409 50.36650
y  34.07928 71.87206
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Grid attributes:
       cellcentre.offset   cellsize cells.dim
x         -31.31151 0.06515996      1254
y          34.11186 0.06515996       580
Data attributes:
Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.11290  0.00000  0.00000  0.06136  0.00000  1.90300

Any thoughts??

Nav
  • 135
  • 2
  • 7

1 Answers1

1

You are getting this message because after reprojection the points that where on a regular grid in lat-lon are no longer a regular grid in another projected system. The reprojection step does not preserve the regularness because for a lot of projections the change between lat-lon and a projected system is not constant. So the distortion, or change, is larger in one area than in another, leading to a non-regular grid. Imagine the squares of the grid in lat-long becoming a non-square shape in a projected system.

The solution that is most simple is just to project the points to EPSG:4326 first, and then perform the interpolation step. I'm not sure which interpolation routine you use, but e.g. gstat does not support interpolation in lat-lon. So, going to a projected system, not lat-lon, before interpolating the data is always a safe bet.

If you really need to reproject a grid, some kind of interpolation needs to take place. You do two things: using a tool like gdalwarp that can perform the coordinate transform and the subsequent warping of the resulting grid. Alternatively, you can interpolate the non-regular grid to a regular grid using e.g. nearest neighbor or any interpolation routine you like.

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
  • Thanks Paul, for your quick reply. The interpolation was carried out on a grid projected in EPSG 3035. Does that mean i should also re project the grid as well to perform the interpolation? Then it will be like coming back to the question again? how do i reproject the grid? – Nav Mar 06 '13 at 21:55
  • It is easiest just to project all the point datasets to the target projection, this is the least hassle. I'll edit some more info into my answer above. – Paul Hiemstra Mar 06 '13 at 22:07
  • Please correct me if i am wrong. The dataset itself is the interpolated values with Lon Lat and corresponding values, which can be converted to a regular grid. If i wish to re-project the interpolated grid, should i use another interpolation routine like idw () in gstat again on the interpolated values. – Nav Mar 06 '13 at 22:29
  • I'd dump the file to a file, e.g. geotiff, and use gdalwarp to change the projection of the grid. – Paul Hiemstra Mar 06 '13 at 22:37
  • is it possible in R, any link with an example. I can convert the dataset to a raster initially and then change the projection to EPSG:4326 using projectRaster and convert it into a geotiff. But how do i use the gdalwarp in R – Nav Mar 06 '13 at 22:40
  • You can call gdalwarp from R using the system command. You could create a function which takes the grid data, dumps it to a temporary file, calls gdalwarp and reads the result. – Paul Hiemstra Mar 06 '13 at 22:53