Skip to content

CRS for the European Biogeographical Regions in R

8 messages · Barbara Szabó, Barry Rowlingson, Roland Kaiser

#
Hi,
i`d be interested if somebody would help me in the CRS for the European
Biogeographical Regions data base.

I have coordinates, already saved as SpatialPolygon in lat-long

sp_poly at bbox

     min      max

x 40.933 59.30397

y 13.633 27.53300
I have to check in which polygons this coordinates belong to. The function
over() should do the job, if both data sets are in the same CRS.


The shape file of the biogeographical regions can be downloaded at
http://www.eea.europa.eu/data-and-maps/data/biogeographical-regions-europe

The CRS is NA when reading it with
library(maptoos)
eea <- readShapeLines("BiogeoRegions2011")
The data lays in the bbox of 943619 489768 8663038 7880872
On the website there is an information on the CRS:
urn:ogc:def:crs:EPSG:6.11:4326

The problem seems to be to specify the correct CRS in eea. First I do not
know the CRS and second if I assume one, I get an error:
proj4string(eea) <- CRS("+proj=longlat +ellps=WGS84")
Error in proj4string. Geographical CRS given to non-conformance data:
8663038.09057 ....

Any ideas which bringing me nearer to a solution are very welcome - thanks
a lot in advance.

Best,
Barbara
#
Hi Barbara,

the data set projection is EPSG:3035 - ETRS89 / ETRS-LAEA.

llibrary(sp)
library(rgdal)

#	the projection information is given in the *.qpj file
#	unfortunately, it is not read from the dataset by the OGR driver
#	we have do define it manually by specifing argument p4s
pg <- readOGR(dsn = "BiogeoRegions2011_shapefile",
	layer = "BiogeoRegions2011", p4s = "+init=epsg:3035")
proj4string(pg)

#	example points
pt <- data.frame(x = c(14,28), y = c(41, 59), id = 1:2)
coordinates(pt) <- ~x+y
proj4string(pt) <- CRS("+init=epsg:4326")
proj4string(pt)

#	project to CRS of pg	
pt <- spTransform(pt, CRS("+init=epsg:3035"))

over(pt, pg)

Cheers, Roli
#
Wow, thank you so, it works!
You saved my day :), otherwise I could spend some more hours on it...

Bests
Barbara


2015-02-07 17:53 GMT+01:00 Roland Kaiser <kardinal.eros at gmail.com>:

  
  
#
On Sat, Feb 7, 2015 at 4:53 PM, Roland Kaiser <kardinal.eros at gmail.com> wrote:

            
I can't see a *.qpj file in that shapefile. There is a *.prj file
which contains projection info.
The *.prj file **is** read when using readOGR:

 > r=readOGR(".","BiogeoRegions2011")
 OGR data source with driver: ESRI Shapefile
 Source: ".", layer: "BiogeoRegions2011"
 with 12 features and 4 fields
 Feature type: wkbPolygon with 2 dimensions
 > proj4string(r)
 [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+ellps=GRS80 +units=m +no_defs"

- that's a full projection specification, rather than an epsg code
though, because that's what's in the file.
No we don't!
Is the .prj file text coding the same projection as epsg:3035?
http://epsg.io/ says 3035 is:

+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

So yes. But there's no need to specify it in readOGR since you have
the .prj file. The maptools functions readShapeSpatial and friends
ignore the .prj file - I don't know why they exist any more!

 > r2=readShapeSpatial("BiogeoRegions2011.shp")
 > proj4string(r2)
[1] NA

Ouch.

Barry
#
Hi!

Sorry for the confusion, I followed the link on the top the page "Note: new version is available!?.
For some reason, this ZIP archive has a *.qpj instead of *.prj file.

Of course, readOGR treats the *prj correctly, if present!

-Roli
#
Dear Roli, dear Barry, dear all,

unfortunately, I cannot find out what is going wrong with my example.
The problem seems to be very basic, but after searching and reading a
lot, I do not come to a solution. It seems that either the
transformation of the WGS84 lat/long to epsg:4326 is somehow wrong or
the biogeographical region CRS is in another reference system than the
.prj file reported. Any help would be wonderful. Thanks in advance.

library(sp)
library(rgdal)

### possible basic problem:
### data in lat/lon -- WGS84, some stations in slovenia:
pt <- data.frame("lat"=c(45.900,46.500,46.067),
"lon"=c(13.633,13.717,14.517), "stid"=c("psi004","psi006","psi003"))
coordinates(pt) <- ~lat+lon
proj4string(pt) <- CRS("+init=epsg:4326")
proj4string(pt)
#[1] "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+towgs84=0,0,0"
## projection seems to be wrong since of negative numbers?
pt at coords

#         lat        lon
#[1,] 8295839 -96345.163
#[2,] 8355036 -58479.199
#[3,] 8286066   1886.133



#################################################
### full example:
pg <- readOGR(dsn = "data/eea",
              layer = "BiogeoRegions2011", p4s = "+init=epsg:3035")
proj4string(pg)

# [1] "+init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000
+y_0=3210000 #+ellps=GRS80 +units=m +no_defs"

## three stations from slovenia:

pt <- data.frame("lat"=c(45.900,46.500,46.067),
"lon"=c(13.633,13.717,14.517), "stid"=c("psi004","psi006","psi003"))

coordinates(pt) <- ~lat+lon
proj4string(pt) <- CRS("+init=epsg:4326")
proj4string(pt)

#[1] "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+towgs84=0,0,0"

##       project to CRS of pg
pt <- spTransform(pt, CRS("+init=epsg:3035"))
res <- over(pt, pg)
res

#  NAME ABBRE code label
#1 <NA>  <NA> <NA>  <NA>
#2 <NA>  <NA> <NA>  <NA>
#3 <NA>  <NA> <NA>  <NA>

## example of a polygon:
summary(pg at polygons[[1]]@Polygons[[1]]@coords)

#       V1                V2
# Min.   :4389214   Min.   :1178199
# 1st Qu.:4396641   1st Qu.:1186794
# Median :4407574   Median :1191444
# Mean   :4405907   Mean   :1193508
# 3rd Qu.:4413052   3rd Qu.:1200195
# Max.   :4421392   Max.   :1207598

pt at coords

#         lat        lon
#[1,] 8295839 -96345.163
#[2,] 8355036 -58479.199
#[3,] 8286066   1886.133



2015-02-08 13:54 GMT+01:00 Roland Kaiser <kardinal.eros at gmail.com>:

  
  
#
Hi Barbara!

the problem is related as to how you specified the formula.
?coordinates says ? in the form of e.g. ~x+y ? (in your case x = lon[gitude], y = lat[itude})

coordiantes(pt) <- lon+lat

is what you want.

?Roli

  
  
#
Hi Roli,

yeah, this was really the problem and now i have what i wanted.
Thank you!

Barbara


2015-02-08 15:26 GMT+01:00 Roland Kaiser <kardinal.eros at gmail.com>: