[R-sig-Geo] CRS for the European Biogeographical Regions in R

[Barbara Szabó]

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!
>
> 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
>
> > Am 08.02.2015 um 12:37 schrieb Barry Rowlingson <
> b.rowlingson at lancaster.ac.uk>:
> >
> > On Sat, Feb 7, 2015 at 4:53 PM, Roland Kaiser <kardinal.eros at gmail.com>
> wrote:
> >
> >> #       the projection information is given in the *.qpj file
> >
> > I can't see a *.qpj file in that shapefile. There is a *.prj file
> > which contains projection info.
> >
> >> #       unfortunately, it is not read from the dataset by the OGR driver
> >
> > 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.
> >
> >> #       we have do define it manually by specifing argument p4s
> >
> > No we don't!
> >
> >> pg <- readOGR(dsn = "BiogeoRegions2011_shapefile",
> >>        layer = "BiogeoRegions2011", p4s = "+init=epsg:3035")
> >
> > 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]