[R-sig-Geo] Help when moving from PROJ4 to PROJ6
Edzer Pebesma
edzer@pebe@m@ @end|ng |rom un|-muen@ter@de
Wed Jul 29 16:01:20 CEST 2020
This seems to have been resolved now by reinstalling sf.
Re: axis order: sf will not swap coordinates by itself, it is only the
way x/y coordinates are interpreted that changes (for EPSG:4326 as
lng/lat by default, but as lat/lng after giving st_axis_order(TRUE)).
You need to do something to get them interpreted as lat/lng.
On 7/28/20 3:34 PM, Roger Bivand wrote:
> On Tue, 28 Jul 2020, Gilberto Camara wrote:
>
>> Dear R-SIG-GEO (esp. Roger and Edzer)
>>
>> I am having problems with “raster”, “rgdal”, and “sf” when moving from
>> PROJ4 to PROJ6. Roger’s explanation on the r-spatial blog
>> ("https://www.r-spatial.org/r/2020/03/17/wkt.html”) and the rgdal blog
>> (http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html
>> <http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html>) are
>> clear on the reasons why there are problems with CRS. As Roger points
>> out, “package maintainers will need to review their use of crs
>> objects”. So very true!
>>
>> I appreciate the tremendous effort of Edzer and Roger for moving from
>> PROJ4 to PROJ6.
>>
>> However, I am failing to understand what is happening in the following
>> example:
>>
>> ================================================
>> library(sf)
>
> Could you please report:
>
> sf_extSoftVersion()
>
> Mine are (on Linux, own PROJ & GDAL installations):
>
> sf_extSoftVersion()
> GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H
> "3.8.1" "3.1.2" "7.1.0" "true" "true"
>
> where there is no change in the coordinates (there should not be a
> change as sf::st_crs() only sets the CRS (but may swap the axes as per
> specification).
>
> I see:
>
>> st_crs(ll_sfc2)
> Coordinate Reference System:
> User input: +proj=longlat +datum=WGS84 +no_defs
> wkt:
> GEOGCRS["unknown",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]],
> ID["EPSG",6326]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8901]],
> CS[ellipsoidal,2],
> AXIS["longitude",east,
> ORDER[1],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> AXIS["latitude",north,
> ORDER[2],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]]]
>
> but
>
>> st_crs(ll_sfc)
> Coordinate Reference System:
> User input: EPSG:4326
> wkt:
> GEOGCRS["WGS 84",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433]],
> CS[ellipsoidal,2],
> AXIS["geodetic latitude (Lat)",north,
> ORDER[1],
> ANGLEUNIT["degree",0.0174532925199433]],
> AXIS["geodetic longitude (Lon)",east,
> ORDER[2],
> ANGLEUNIT["degree",0.0174532925199433]],
> USAGE[
> SCOPE["unknown"],
> AREA["World"],
> BBOX[-90,-180,90,180]],
> ID["EPSG",4326]]
>
> I can't see from
> https://www.r-spatial.org/r/2020/03/17/wkt.html#axis-order how to
> instantiate a crs object with non-authorized order.
>
> In sp/rgdal:
>
>> sp <- SpatialPoints(cbind(longitude, latitude),
> proj4string=CRS("+proj=longlat +datum=WGS84"))
>> sp
> SpatialPoints:
> longitude latitude
> [1,] -55.62277 -11.75932
> Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
> +no_defs
>> cat(wkt(sp), "\n")
> GEOGCRS["unknown",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]],
> ID["EPSG",6326]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8901]],
> CS[ellipsoidal,2],
> AXIS["longitude",east,
> ORDER[1],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> AXIS["latitude",north,
> ORDER[2],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]]]
>
>> sp2 <- SpatialPoints(cbind(longitude, latitude),
> proj4string=CRS(SRS_string="EPSG:4326"))
>> sp2
> SpatialPoints:
> longitude latitude
> [1,] -55.62277 -11.75932
> Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
> +no_defs
>> cat(wkt(sp2), "\n")
> GEOGCRS["WGS 84",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]],
> ID["EPSG",6326]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8901]],
> CS[ellipsoidal,2],
> AXIS["longitude",east,
> ORDER[1],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> AXIS["latitude",north,
> ORDER[2],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> USAGE[
> SCOPE["unknown"],
> AREA["World"],
> BBOX[-90,-180,90,180]]]
>
> where sp/rgdal enforces legacy/GIS/visualization axis order. Coercing
> from the EPSG:4326 sp/rgdal object to sfc gives the intended order:
>
>> ll_sfc3 <- st_as_sfc(sp2)
>> ll_sfc3
> Geometry set for 1 feature
> geometry type: POINT
> dimension: XY
> bbox: xmin: -55.62277 ymin: -11.75932 xmax: -55.62277 ymax:
> -11.75932
> geographic CRS: WGS 84
> POINT (-55.62277 -11.75932)
>> st_crs(ll_sfc3)
> Coordinate Reference System:
> User input: WGS 84
> wkt:
> GEOGCRS["WGS 84",
> DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]],
> ID["EPSG",6326]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8901]],
> CS[ellipsoidal,2],
> AXIS["longitude",east,
> ORDER[1],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> AXIS["latitude",north,
> ORDER[2],
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> USAGE[
> SCOPE["unknown"],
> AREA["World"],
> BBOX[-90,-180,90,180]]]
>
> @Edzer - should we raise an sf issue, as st_crs<- appears not to respect
> st_axis_order()?
>
> I'm not sure that this clarifies enough, but it's a start ...
>
> Roger
>
>>> longitude <- -55.62277
>>> latitude <- -11.75932
>>
>>> st_point <- sf::st_point(c(longitude, latitude))
>>> ll_sfc <- sf::st_sfc(st_point, crs = "EPSG:4326")
>>
>>> ll_sfc[[1]]
>> POINT (-55.55527 -11.51782)
>>
>>> ll_sfc2 <- sf::st_sfc(st_point, crs = "+proj=longlat +datum=WGS84
>>> +no_defs”)
>>
>>> ll_sfc2[[1]]
>> POINT (-55.62277 -11.75932)
>> ==================================================
>>
>>> sessionInfo()
>> R version 4.0.1 (2020-06-06)
>> Platform: x86_64-apple-darwin17.0 (64-bit)
>> Running under: macOS Catalina 10.15.5
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] sf_0.9-5 sits_0.9.5.1 raster_3.3-13 rgdal_1.5-12 sp_1.4-2
>> =====================================================================
>>
>> Why is using the EPSG code resulting in a different behaviour that
>> using the old PROJ4 text?
>>
>> Many thanks
>> Gilberto
>>
>> ===========================
>> Prof Dr Gilberto Camara
>> Secretariat Director
>> GEO - Group on Earth Observations
>> 7 bis, Avenue de La Paix
>> CH-1211 Geneva - Switzerland
>> Tel: +41227308480
>> www.earthobservations.org
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo using r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48149 Muenster, Germany
Phone: +49 251 8333081
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 3110 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20200729/5bb3203e/attachment.bin>
More information about the R-sig-Geo
mailing list