Help when moving from PROJ4 to PROJ6

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Help when moving from PROJ4 to PROJ6

Gilberto Camara
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)
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Help when moving from PROJ4 to PROJ6

Roger Bivand
Administrator
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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Help when moving from PROJ4 to PROJ6

Gilberto Camara
Dear Roger

Please see my output:

> sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H
       "3.8.1"        "3.1.1"        "6.3.1"         "true"         “true"

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

> On 28 Jul 2020, at 15:34, Roger Bivand <[hidden email]> wrote:
>
> sf_extSoftVersion()


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Help when moving from PROJ4 to PROJ6

Gilberto Camara
Dear Roger

I have installed the latest version of project, using “brew” as recommended by the PROJ website. I have also installed the latest versions of “sf” and “rgdal” as recommended by the “sf” GitHub.

install.packages("rgdal", configure.args = c("--with-proj-lib=/usr/local/lib/", "--with-proj-include=/usr/local/include/“))

However, I still seeing the same environment for “sf” (proj version 6.3.1). Is there an additional command I need to do?

Thanks
Gilberto


> On 28 Jul 2020, at 15:36, Gilberto Camara <[hidden email]> wrote:
>
> Dear Roger
>
> Please see my output:
>
>> sf_extSoftVersion()
>          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H
>       "3.8.1"        "3.1.1"        "6.3.1"         "true"         “true"
>
> 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
>
>> On 28 Jul 2020, at 15:34, Roger Bivand <[hidden email]> wrote:
>>
>> sf_extSoftVersion()
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Help when moving from PROJ4 to PROJ6

edzer
In reply to this post by Roger Bivand
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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48149 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (3K) Download Attachment