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.biv...@nhh.no>
Enviado: quarta-feira, 16 de junho de 2021 12:01
Para: Pietro Andre Telatin Paschoalino <pietro_tel...@hotmail.com>
Cc: r-sig-geo@r-project.org <r-sig-geo@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.biv...@nhh.no>
> Enviado: quarta-feira, 16 de junho de 2021 08:38
> Para: Pietro Andre Telatin Paschoalino <pietro_tel...@hotmail.com>
> Cc: r-sig-geo@r-project.org <r-sig-geo@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.biv...@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.biv...@nhh.no>
>> > Enviado: ter�a-feira, 15 de junho de 2021 15:18
>> > Para: Pietro Andre Telatin Paschoalino
>> <pietro_tel...@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.biv...@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.biv...@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.biv...@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.biv...@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
