Hi List,

I am using the rgdal package to project lat/lon co-ordinates to utm.  I
then use an interpolation function to generate a new set of utm
co-ordinates, which I need to re-project as lat/lon.  I've been testing
spTransform() and project() from rgdal for re-projection, using original
co-ordinates.

I encounter errors when projecting using the correct zone
http://spatialreference.org/ref/epsg/32630/  (for example: 51.73, -5.47),
but these errors are much diminished if I use a different zone:

ERRORS, CORRECT ZONE:
> g151_sp <- SpatialPoints(g151_latlon, proj4string=CRS("+proj=longlat
+zone=30N +datum=WGS84 +ellps=WGS84"))
> g151_tr_utm <- spTransform(g151_sp, CRS("+proj=utm +zone=30N +datum=WGS84
+ellps=WGS84 +units=m +no_defs"))
> g151_utm <- coordinates(g151_tr_utm)
> head(g151_utm)
     coords.x1 coords.x2
[1,]   7727626  -1040054
[2,]   7727999  -1038178
[3,]   7728274  -1037516
[4,]   7728477  -1037170
[5,]   7728336  -1037251
[6,]   7728312  -1037115
> g151_tr_ll <- SpatialPoints(g151_tr_utm, proj4string=CRS("+proj=utm
+zone=30N +datum=WGS84 +ellps=WGS84 +units=m +no_defs"))
> g151_ll <- spTransform(g151_tr_utm, CRS("+proj=longlat +zone=30N
+datum=WGS84 +ellps=WGS84"))
> head(g151_ll)
SpatialPoints:
     coords.x1 coords.x2
[1,]  50.92425 -5.779419
[2,]  50.92546 -5.768819
[3,]  50.92652 -5.765044
[4,]  50.92734 -5.763052
[5,]  50.92674 -5.763544
[6,]  50.92660 -5.762785
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0
> head(g151_latlon)
         [,1]     [,2]
[1,] 51.73138 -5.47851
[2,] 51.73208 -5.46843
[3,] 51.73309 -5.46478
[4,] 51.73393 -5.46283
[5,] 51.73324 -5.46335
[6,] 51.73302 -5.46265
LESS ERROR, WRONG ZONE:
> g151_sp <- SpatialPoints(g151_latlon, proj4string=CRS("+proj=longlat
+zone=36N +datum=WGS84 +ellps=WGS84"))
> g151_tr_utm <- spTransform(g151_sp, CRS("+proj=utm +zone=36N +datum=WGS84
+ellps=WGS84 +units=m +no_defs"))
> g151_utm <- coordinates(g151_tr_utm)
> head(g151_utm)
     coords.x1 coords.x2
[1,]   2612390 -639448.1
[2,]   2612509 -638275.0
[3,]   2612641 -637853.0
[4,]   2612747 -637628.7
[5,]   2612664 -637686.8
[6,]   2612641 -637604.3
> g151_tr_ll <- SpatialPoints(g151_tr_utm, proj4string=CRS("+proj=utm
+zone=36N +datum=WGS84 +ellps=WGS84 +units=m +no_defs"))
> g151_ll <- spTransform(g151_tr_utm, CRS("+proj=longlat +zone=36N
+datum=WGS84 +ellps=WGS84"))
> head(g151_ll)
SpatialPoints:
     coords.x1 coords.x2
[1,]  51.73135 -5.478509
[2,]  51.73205 -5.468429
[3,]  51.73306 -5.464779
[4,]  51.73390 -5.462829
[5,]  51.73321 -5.463349
[6,]  51.73299 -5.462649
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0
> head(g151_latlon)
         [,1]     [,2]
[1,] 51.73138 -5.47851
[2,] 51.73208 -5.46843
[3,] 51.73309 -5.46478
[4,] 51.73393 -5.46283
[5,] 51.73324 -5.46335
[6,] 51.73302 -5.46265

Exactly the same happens if I use project:
g151_proj <- project(g151_latlon, "+proj=utm +datum=WGS84 +zone=36N
+ellps=WGS84 +units=m", inv=F)
g151_reproj <- project(g151_proj, "+proj=utm +datum=WGS84 +zone=36N
+ellps=WGS84 +units=m", inv=T)

Any guidance much appreciated,


Kylie Scales
Postgraduate Researcher
Plymouth Marine Laboratory
+44 (0)1752 633165

	[[alternative HTML version deleted]]

