[R-sig-Geo] How to translate coordinates from cylindrical to a Lambert Azimuthal Equal Area Projection

Barry Rowlingson b.rowlingson at lancaster.ac.uk
Tue Apr 23 18:41:31 CEST 2013


Okay, maybe all you want to do is work out where your raster projects to?

Here's a raster in epsg:3035 - you have something like this?

 > re
class       : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 7785.786, 11293.11  (x, y)
extent      : 3577016, 3654874, 3714447, 3827378  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10
+x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
data source : in memory
names       : layer
values      : 0.001494624, 0.9948669  (min, max)

Now, to see where that is in lat-long, simply do:

 > extent(projectExtent(re,CRS("+init=epsg:4326")))
class       : Extent
xmin        : -2.290917
xmax        : -0.7665775
ymin        : 56.07447
ymax        : 57.08626

(But as I keep banging on, your grid in epsg:3035 won't line up with
lat-long lines in that box)

The other thing you can do is use spTransform to work out where the
four corners of your raster transform to.

This corners function gets the four corner coordinates of a raster:

corners <- function(r){
  e = extent(r)
  x=c(xmin(e),xmin(e),xmax(e),xmax(e))
  y=c(ymin(e),ymax(e),ymin(e),ymax(e))
  SpatialPoints(data.frame(x=x,y=y),proj4string=CRS(projection(r)))
}

and returns a SpatialPoints object in the original coordinate system:

 > corners(re)
SpatialPoints:
           x       y
[1,] 3577016 3714447
[2,] 3577016 3827378
[3,] 3654874 3714447
[4,] 3654874 3827378
Coordinate Reference System (CRS) arguments: +init=epsg:3035 +proj=laea
+lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m
+no_defs

Let's see where these transform to:

 > spTransform(corners(re),CRS("+init=epsg:4326"))
SpatialPoints:
             x        y
[1,] -1.972778 55.97462
[2,] -2.290917 56.97278
[3,] -0.738697 56.08479
[4,] -1.025412 57.08626
Coordinate Reference System (CRS) arguments: +init=epsg:4326
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

Useful?

Barry



More information about the R-sig-Geo mailing list