Tuesday 17 February 2015

r - What is the proj4string of the Lambert azimuthal equal-area projection?


Shapefiles of European data are available on the Eurostat website. I have loaded the NUTS_2013_60M_SH Shapefile in R with the command;


europeRG <- readOGR(dsn = patheurope, layer = "NUTS_RG_60M_2013", stringsAsFactors = FALSE)

I have managed to transform the map of Europe to the web mercator projection system used by OpenStreetMap and Google maps (EPSG:3857).


europewmercator <- spTransform(europeRG, CRS("+init=epsg:3857"))
plot(europewmercator)


But a 2001 JRC and EuroGeographics document entitled Map projections for Europe recommends using the Lambert Azimuthal Equal Area projection:



[...] Based on the above, it is documented that the projection to be used for statistical mapping in the EU, is the Lambert Azimuthal Equal Area projection. The projection center should be positioned at 9° E and 53° N. [...]



Wikipedia page of the Lambert azimuthal equal-area projection. This site says the epsg code of the "Lambert Azimuthal Equal Area" projection is 9820. Unfortunately EPSG code 9820 is not available on https://epsg.io/?q=9820 or http://spatialreference.org but it is available on http://www.epsg-registry.org/ The issue is that the R function sptransform() doesn't recognise epsg:9820.


plot(europeRG)
# Transform To Lambert Azimuthal Equal Area projection
europeLAEA <- spTransform(europeRG, CRS("+init=epsg:9820"))


Error in spTransform(europeRG, CRS("+init=epsg:9820")) :
error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': Error in CRS("+init=epsg:9820") : no options found in 'init' file

What proj4string could I use to perform the conversion?



Answer



9820 is the EPSG code for the laea projection definition. You can find it at https://www.epsg-registry.org/ under retrieve by code. This code numer has nothing to do with the commonly used EPSG codes like 4326 and 3857. These include projection parameters specified for a certain country or region.


For laea centered on Europe, the code is EPSG:3035, with these parameters:


+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

As you can see, the center is at 10°E and 52°N, which is not exactly what your source recommends. If you want those, you can build a custom CRS with parameters:



+proj=laea +lat_0=53 +lon_0=9 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The x_0 and y_0 have no effect on the projection, they just shift the coordinates to avoid negative coordinates.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...