Projection systems
From Intamap
Contents |
Projections
Projection and reprojection can be done using ArcGIS (yes, obviously) but also using the open source library proj4. Luckily, many projects linked to proj4, such as R (package rgdal, function spTransform) and PostGIS
We distinguish geographical coordinates (spherical: longitude/latitude in decimal degrees or DMS = DegreesMinutesSeconds, 53,42'33") and projected coordinates. Doing spatial statistics in geographical coordinates is hindered by the fact that
- we need special formula's for distances (e.g. the great spherical distance, shortest path measured along the sphere) and
- that our data are not 2D, implying we need to make sure that covariance functions are positive definite on the sphere
PROJ4 is a library for cartographic (re)projections. It is interfaced from many open source projects, including R (package rgdal) and PostGIS. It includes the EPSG set of projections, two of which are relevant here:
Geographic coordinates
+init=epsg:4326 which is short for projection ID WGS 84: full parameters:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
(meaning: geographical coordinates, i.e. longitude (x) latitude (y) in decimal degrees, using the WGS84 spheroid as its base; this spheroid seems to be the default if we don't know better)
INSPIRE projection
+init=epsg:3035 which is short for projection ID ETRS-LAEA (Lambert Equal Area, ETRS89 geodetic datum, GRS80 ellipsoid); full parameters:
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
INSPIRE-compliant reprojection in R
Sounds good, let's try. After insalling R with packages sp and rgdal, try
library(rgdal)
tst = SpatialPoints(data.frame(x=c(5,5),y=c(50,60)))
proj4string(tst) = CRS("+init=epsg:4326")
tst
spTransform(tst, CRS("+init=epsg:3035"))
and the output is
> library(rgdal)
Loading required package: sp
Geospatial Data Abstraction Library extensions to R successfully loaded
> tst = SpatialPoints(data.frame(x=c(5,5),y=c(50,60)))
> proj4string(tst) = CRS("+init=epsg:4326")
> tst
SpatialPoints:
x y
[1,] 5 50
[2,] 5 60
Coordinate Reference System (CRS) arguments: +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
> spTransform(tst, CRS("+init=epsg:3035"))
SpatialPoints:
x y
[1,] 3962799 2999719
[2,] 4041548 4109792
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
which reproduces Table 2 example (p 26) in Annoni's report below, when rounded to metres.
The following
library(maps)
library(mapdata)
library(maptools)
wrld = map("world", interior = FALSE, plot = FALSE)
wrld = pruneMap(wrld)
wrld.sp = map2SpatialLines(wrld, proj4string = CRS("+init=epsg:4326"))
INSP.prj = CRS("+init=epsg:3035")
wrld_grd = gridlines(wrld.sp, easts=seq(-180,160,10), norths=seq(-50,90,10), ndiscr=100)
library(rgdal)
wrld.sp = spTransform(wrld.sp, INSP.prj)
wrld_grd = spTransform(wrld_grd, INSP.prj)
plot(wrld_grd, xlim= c(0, 1e7), ylim = c(0, 1e7), col = 'grey')
lines(wrld.sp)
box()
title("Figure 1 (p13) Annoni report")
produces
which bears a striking resemblance to figure 1 in the Annoni report, see below.
Grids
The report that describes definition is European Reference Grids Workshop, edited by A. Annoni, European Commission Joint Research Centre, EUR 21494 EN. You can download shapefiles from the reference grids (1km, 10km, 100km) from e.g. the eea or from jrc (found by googling on "european reference grids").
The grids are defined on powers of 10, unit m: 1m, 10m, 100m, 1km, 10km, 100km, 1000km, hierarchically (meaning they have a common origin). Also, up to two quad-tree subdivisions are allowed, e.g. 500m and 250m, or 50km and 25 km, but not 12.5km (use 10km instead). The report then explains how for each grid cell a unique code can be derived from the coordinates of the south-western corner of the grid cell.
Coming: how to generate such a grid in R.
