[R-sig-Geo] 3 Questions on CRS change

Pietro Andre Telatin Paschoalino p|etro_te|@to @end|ng |rom hotm@||@com
Wed Jun 16 17:38:02 CEST 2021


Thank you Roger. I think that the difference between the two will not be huge but I'll test with this s2 package.

Only a doubt, how to transform a shapefile opened by readOGR or st_read to a as_s2_geography (I can't extract the centroids in those formats in s2 package),
and when I try to transform I got a message:

Error in UseMethod("as_s2_geography") :
 no applicable method for 'as_s2_geography' applied to an object of class "c('sf', 'data.frame')"

Thanks again.

Pietro Andre Telatin Paschoalino
Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.

________________________________
De: Roger Bivand <Roger.Bivand using nhh.no>
Enviado: quarta-feira, 16 de junho de 2021 12:01
Para: Pietro Andre Telatin Paschoalino <pietro_telato using hotmail.com>
Cc: r-sig-geo using r-project.org <r-sig-geo using r-project.org>
Assunto: RE: [R-sig-Geo] 3 Questions on CRS change

On Wed, 16 Jun 2021, Pietro Andre Telatin Paschoalino wrote:

> Hello Roger,
>
> After updating the packages issue 1 was solved.
>

Many data sets in the wild declare EPSG:4326 but use eastings:northings
contrary to the authority. sp/rgdal/raster use eastings/northings and
rgdal changes the recorded order in the crs/CRS object (used by
sp/raster). OGC:CRS84 is eastings/northings by definition. Caution is
advised. Further, from PROJ 8 and EPSG version >= 10, WGS84 is a datum
ensemble:

> sp::CRS("OGC:CRS84")
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["WGS 84 (CRS84)",
     ENSEMBLE["World Geodetic System 1984 ensemble",
         MEMBER["World Geodetic System 1984 (Transit)"],
         MEMBER["World Geodetic System 1984 (G730)"],
         MEMBER["World Geodetic System 1984 (G873)"],
         MEMBER["World Geodetic System 1984 (G1150)"],
         MEMBER["World Geodetic System 1984 (G1674)"],
         MEMBER["World Geodetic System 1984 (G1762)"],
         ELLIPSOID["WGS 84",6378137,298.257223563,
             LENGTHUNIT["metre",1]],
         ENSEMBLEACCURACY[2.0]],
     PRIMEM["Greenwich",0,
         ANGLEUNIT["degree",0.0174532925199433]],
     CS[ellipsoidal,2],
         AXIS["geodetic longitude (Lon)",east,
             ORDER[1],
             ANGLEUNIT["degree",0.0174532925199433]],
         AXIS["geodetic latitude (Lat)",north,
             ORDER[2],
             ANGLEUNIT["degree",0.0174532925199433]],
     USAGE[
         SCOPE["Not known."],
         AREA["World."],
         BBOX[-90,-180,90,180]],
     ID["OGC","CRS84"]]
> rgdal::rgdal_extSoftVersion()
           GDAL GDAL_with_GEOS           PROJ             sp           EPSG
        "3.3.0"         "TRUE"        "8.0.1"        "1.4-6"      "v10.018"

where it is expected that current use will assume (G1762) rather than
earlier versions. This impacts especially parts of the world with
fast-moving tectonic plates. We don't yet know how it will play out, nor
how to choose non-default members of the ensemble.

> 2- I keep my doubt if the order axis between the shape and raster (or 2
>    shapes) in cat(wkt(), "\n", before the spTransform is a problem. But
>    I don't think so, I tested it with another shape without spTransform
>    and the results were similar. I believe spTransform also changes the
>    axis order correctly. But I don't know if it's a specific case or if
>    it won't always be a problem to apply spTransform in these cases,
>    would you know about this?
>
> 3- In this case, as I don't know what you say by unprojected in my case,
>    you believe it's a problem to apply gCentroid in my example:

Unprojected is GEOGCRS, projected is PROJCRS (plus some other variants).
For GEOS to work, you need planar, projected coordinates where Euclidean
distances work, rather than distances on the surface of a sphere. The s2
package provides s2::s2_centroid() on the surface of a sphere.

Roger

>
> 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]]]]
>
> Thank you.
>
> Pietro Andre Telatin Paschoalino
> Doutorando em Ciências Econômicas da Universidade Estadual de Maringá - PCE.
>
> ________________________________
> De: Roger Bivand <Roger.Bivand using nhh.no>
> Enviado: quarta-feira, 16 de junho de 2021 08:38
> Para: Pietro Andre Telatin Paschoalino <pietro_telato using hotmail.com>
> Cc: r-sig-geo using r-project.org <r-sig-geo using r-project.org>
> Assunto: Re: [R-sig-Geo] 3 Questions on CRS change
>
> On Wed, 16 Jun 2021, Pietro Andre Telatin Paschoalino wrote:
>
>> Hello Roger,
>> Thank you. I'll update the packages and try again.
>>
>> And about my doubts 2 and 3, you know If there is any problem?
>
> Try them after updating as they may differ. 2 - nobody knows,
> sp/rgdal/raster try to impose GIS/visualization order. 3 - try to avoid
> using rgeos with unprojected geometries, it only provides accurate
> centroids for planar geometries.
>
> Roger
>
>>
>> Thanks again.
>>
>> Em 16 de jun de 2021 04:18, Roger Bivand <Roger.Bivand using nhh.no> escreveu:
>>
>>       On Tue, 15 Jun 2021, Pietro Andre Telatin Paschoalino wrote:
>>
>>      > 1-  I tried to compare the CRS OGC with my raster (stack)
>>       and I couldn't:
>>      >
>>      > My raster data are CRU DATA TS 4.04 (EAST ANGLIA)
>>
>>       This is not a minimal reproducible example, perhaps
>>
>>       https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04/
>>
>>       is the source, but no specific (very small) file is indicated.
>>       Using the
>>       224MB cru_ts4.04.2011.2019.pre.dat.nc on an updated system
>>       (rgdal, raster
>>       and sp as CRAN, Linux, PROJ 8.0.1, GDAL 3.3.0) there is no
>>       problem.
>>
>>       You have rgdal 1.5-10 (CRAN has 1.5-23), sp 1.4-2 (CRAN has
>>       1.4-5) and
>>       raster 3.3-13 (CRAN has 3.4-10), which I will not re-install
>>       because they
>>       are not the problem, and the current versions are available
>>       for R 3.6 on
>>       Windows. All the debugging needed is in the output of
>>       sessionInfo().
>>
>>       Please report back after updating at least sp, rgdal and
>>       raster (your
>>       raster version is also pretty outdated).
>>
>>       Roger Bivand
>>
>>      >
>>      > rgdal::compare_CRS(CRS("OGC:CRS84"), slot(rasbrick, "crs"))
>>      > Error in CRS("OGC:CRS84") :
>>      > PROJ4 argument-value pairs must begin with +: OGC:CRS84
>>      >
>>      > I also tried it like this:
>>      >
>>      > rgdal::compare_CRS(CRS("+OGC:CRS84"), slot(rasbrick, "crs"))
>>      > Error in showSRID(uprojargs, format = "PROJ", multiline =
>>       "NO") :
>>      > Can't parse PROJ.4-style parameter string
>>      > +OGC:CRS84
>>      >
>>      > could you answer me where am I wrong?
>>      >
>>      > When I compare a shapefile with my raster (stack) the
>>       function runs normally.
>>      >
>>      > 2- using:
>>      >
>>      > cat(wkt(shapebefore), "\n")
>>      >
>>      > cat(wkt(rasbrick), "\n")
>>      >
>>      > I saw that the order of the axis were different, but after
>>       the transformation they were the same.
>>      >
>>      > Any problem making shapeafter<- spTransform(shapebefore,
>>       crs(rasbrick))
>>      >
>>      > when are the axis different before the transformation?
>>      >
>>      > Ex:
>>      >
>>      > cat(wkt(shapebefore), "\n")
>>      > 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["latitude",north,
>>      >            ORDER[1],
>>      >            ANGLEUNIT["degree",0.0174532925199433]],
>>      >        AXIS["longitude",east,
>>      >            ORDER[2],
>>      >            ANGLEUNIT["degree",0.0174532925199433]],
>>      >    ID["EPSG",4326]]
>>      >
>>      >> cat(wkt(rasbrick), "\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]]]]
>>      >
>>      > The order of the axis (and other parameters) of my new shape
>>       (shapeafter) is the same as the raster (rasbrick) after
>>       transformation.
>>      >
>>      > 3 -  when I use:
>>      >
>>      > centroids = gCentroid(shapebefore, byid=TRUE)
>>      >
>>      > and
>>      >
>>      > centroids = gCentroid(shapeafter, byid=TRUE)
>>      >
>>      > why are the coordinates x and y values equal (and the values
>>       are exactly the same)?
>>      >
>>      > (It seems to me that both are longitude=x and latitude=y)
>>       but since the order is switched in cat(wkt(), shouldn't it be
>>       latitude=x and longitude=y in the shapebefore?
>>      >
>>      > if I use proj4string()
>>      > both my shapebefore and my shapeafter are like:
>>      >
>>      > "+proj=longlat +datum=WGS84 +no_defs"
>>      >
>>      >> sessionInfo()
>>      > R version 3.6.2 (2019-12-12)
>>      > Platform: x86_64-w64-mingw32/x64 (64-bit)
>>      > Running under: Windows 10 x64 (build 19042)
>>      >
>>      > Matrix products: default
>>      >
>>      > locale:
>>      > [1] LC_COLLATE=English_United States.1252
>>       LC_CTYPE=English_United States.1252
>>      > [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
>>      > [5] LC_TIME=English_United States.1252
>>      >
>>      > attached base packages:
>>      > [1] stats     graphics  grDevices utils     datasets
>>       methods   base
>>      >
>>      > other attached packages:
>>      > [1] spdep_1.1-5    sf_0.9-5       spData_0.3.8
>>       maptools_1.0-1 xlsx_0.6.3     rgeos_0.5-3    raster_3.3-13
>>      > [8] rgdal_1.5-10   sp_1.4-2
>>      >
>>      > loaded via a namespace (and not attached):
>>      > [1] Rcpp_1.0.5         compiler_3.6.2     pillar_1.4.4
>>       LearnBayes_2.15.1  class_7.3-15
>>      > [6] tools_3.6.2        boot_1.3-23        nlme_3.1-142
>>       lifecycle_0.2.0    tibble_3.0.1
>>      > [11] lattice_0.20-41    pkgconfig_2.0.3
>>       rlang_0.4.7        Matrix_1.2-18      DBI_1.1.0
>>      > [16] rstudioapi_0.11    expm_0.999-5
>>       coda_0.19-3        rJava_0.9-12       e1071_1.7-3
>>      > [21] dplyr_1.0.0        gtools_3.8.2
>>       generics_0.0.2     xlsxjars_0.6.1     vctrs_0.3.1
>>      > [26] classInt_0.4-3     grid_3.6.2
>>       tidyselect_1.1.0   glue_1.4.1         spDataLarge_0.4.1
>>      > [31] R6_2.4.1           foreign_0.8-72
>>       gdata_2.18.0       deldir_0.1-28      purrr_0.3.4
>>      > [36] magrittr_1.5       gmodels_2.18.1
>>       splines_3.6.2      MASS_7.3-51.4      codetools_0.2-16
>>      > [41] ellipsis_0.3.1     units_0.6-7
>>       KernSmooth_2.23-16 crayon_1.3.4
>>      >
>>      > Thanks.
>>      >
>>      > Pietro Andre Telatin Paschoalino
>>      > Doutorando em Ci�ncias Econ�micas da Universidade Estadual
>>       de Maring� - PCE.
>>      >
>>      > ________________________________
>>      > De: Roger Bivand <Roger.Bivand using nhh.no>
>>      > Enviado: ter�a-feira, 15 de junho de 2021 15:18
>>      > Para: Pietro Andre Telatin Paschoalino
>>       <Pietro_telato using hotmail.com>
>>      > Assunto: Re: Doubt in function compare_CRS
>>      >
>>      > Please do not email me directly. For use questions, please
>>       subscribe to
>>      > R-sig-geo
>>      >
>>      > https://stat.ethz.ch/mailman/listinfo/R-SIG-Geo
>>      >
>>      > and post there. Your problem cannot be reproduced by me (I
>>       do not know
>>      > what your raster is) - when you post provide a minimal
>>       self-contained
>>      > example. Do also provide the output of sessionInfo(), as you
>>       may not be
>>      > using updated packages.
>>      >
>>      > Roger Bivand
>>      >
>>      >
>>      > On Tue, 15 Jun 2021, Pietro Andre Telatin Paschoalino wrote:
>>      >
>>      >> Hello Bivand,
>>      >>
>>      >> How you create the function compare_CRS and it has already
>>       exemplified
>>      >> several issues for me. I would like to know if you can help
>>       me with a
>>      >> simple issue
>>      >>
>>      >> I tried to compare the CRS OGC with my raster (stack) and I
>>       couldn't:
>>      >>
>>      >> rgdal::compare_CRS(CRS("OGC:CRS84"), slot(rasbrick, "crs"))
>>      >> Error in CRS("OGC:CRS84") :
>>      >> PROJ4 argument-value pairs must begin with +: OGC:CRS84
>>      >>
>>      >> I also tried it like this:
>>      >>
>>      >> rgdal::compare_CRS(CRS("+OGC:CRS84"), slot(rasbrick,
>>       "crs"))
>>      >> Error in showSRID(uprojargs, format = "PROJ", multiline =
>>       "NO") :
>>      >> Can't parse PROJ.4-style parameter string
>>      >> +OGC:CRS84
>>      >>
>>      >> could you answer me where am I wrong?
>>      >>
>>      >> When I compare a shapefile with my raster (stack) the
>>       function runs normally, for example:
>>      >>
>>      >> compare_CRS(slot(shape1, "proj4string"), slot(rasbrick,
>>       "crs"))
>>      >>                      strict                   equivalent
>>       equivalent_except_axis_order
>>      >>                       FALSE
>>       TRUE                         TRUE
>>      >>
>>      >> Thank you.
>>      >>
>>      >> Pietro Andre Telatin Paschoalino
>>      >> Doutorando em Ci�ncias Econ�micas da Universidade
>>       Estadual de Maring� - PCE.
>>      >>
>>      >>[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-an
>> imated-no-repeat-v1.gif]<https://www.avast.com/sig-email?utm_medium=emai
>> l&utm_source=link&utm_campaign=sig-email&utm_content=webmail&utm_term=ic
>>       on>    Virus-free.www.avast.com<https://www.avast.com/sig-email?utm_medium=email&utm_source
>>       =link&utm_campaign=sig-email&utm_content=webmail&utm_term=link>
>>      >>
>>      >
>>      > --
>>      > Roger Bivand
>>      > Emeritus Professor
>>      > Department of Economics, Norwegian School of Economics,
>>      > Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
>>      > e-mail: Roger.Bivand using nhh.no
>>      > https://orcid.org/0000-0003-2392-6140
>>      > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>      >
>>      > [[alternative HTML version deleted]]
>>      >
>>      >
>>
>>       --
>>       Roger Bivand
>>       Emeritus Professor
>>       Department of Economics, Norwegian School of Economics,
>>       Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
>>       e-mail: Roger.Bivand using nhh.no
>>       https://orcid.org/0000-0003-2392-6140
>>       https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>>
>>
>>
>
> --
> Roger Bivand
> Emeritus Professor
> Department of Economics, Norwegian School of Economics,
> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
> e-mail: Roger.Bivand using nhh.no
> https://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>

--
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list