[R-sig-Geo] 3 Questions on CRS change
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Wed Jun 16 17:44:51 CEST 2021
On Wed, 16 Jun 2021, Pietro Andre Telatin Paschoalino wrote:
> 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:
With sf 1.0-0 and sf::sf_use_s2() TRUE, sf::st_centroid() should use s2
and convert the objects. Alternatively, use sf::st_as_s2().
Roger
>
> 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
>
--
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
More information about the R-sig-Geo
mailing list