Message-ID: <a5f2defe-8fe0-f4dc-4f8d-f19595adf715@uni-muenster.de>
Date: 2020-07-29T14:01:20Z
From: Edzer Pebesma
Subject: Help when moving from PROJ4 to PROJ6
In-Reply-To: <alpine.LFD.2.23.451.2007281517550.520703@reclus.nhh.no>
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 at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at 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>