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

[Barbara Szabó]

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>:

> 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
> #[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]]