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

[Roland Kaiser]

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


> Am 08.02.2015 um 15:17 schrieb Barbara Szabó <barbara.szabo.elte at gmail.com>:
> 
> 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 <tel: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 <tel: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 <mailto: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 <mailto:b.rowlingson at lancaster.ac.uk>>:
> >
> > On Sat, Feb 7, 2015 at 4:53 PM, Roland Kaiser <kardinal.eros at gmail.com <mailto: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/ <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 <mailto:R-sig-Geo at r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> 
> 
> 
> 


	[[alternative HTML version deleted]]