[R-sig-Geo] Simple doubt - CRS

[Roger Bivand]

On Wed, 24 Mar 2021, Pietro Andre Telatin Paschoalino wrote:

> Hello Roger, thank you for your explanation.
>
> So if my str(crs(rasbrick)) are equal than cat(wkt(shape2), "\n") (here 
> you put shape but I think that is shape2 with the new CRS) so I'm okay, 
> right? There is some function to check automatically if the two are 
> equal?

No, str(crs(rasbrick)) only shows the revealed contents of the CRS object. 
The WKT2-2019 representation has to be hidden in comment(crs(rasbrick)), 
which is *not* shown by str(). It had to be hidden in this way because CRS 
is a formal class, so its definition is fixed. If the formal definition of 
the CRS class was changed, many packages and workflows would break, so the 
WKT2-2019 representation was added, but can most easily be seen using the 
wkt() method.

>
> I'm only don't understand very well what I need to know about the: 
> "get_last_coordOp()" I did that and in return I had:
>
> "+proj=noop +ellps=GRS80"
>

This means that shape and rasbrick had the same CRS, so the coordinate 
operation was a noop, no operation. If they had differed, it would show 
the steps taken. rgdal::list_coordOps() can show what spTransform can 
choose, equivalently in sf 0.9-8 it is called sf::sf_proj_pipelines(). 
https://github.com/r-spatial/sf/issues/1634 is an example of its use. In a 
different case in Canada, not using the PROJ CDN for on-demand 
transformation grid downloading led to 20m errors, so keepin one's eye on 
the instantiable pipelines is potentially important.

Roger

> 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, 24 de março de 2021 11:11
> 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] Simple doubt - CRS
>
> On Wed, 24 Mar 2021, Pietro Andre Telatin Paschoalino wrote:
>
>> Hello everyone,
>>
>> I'm extracting a precipitation data of a raster to shapefile and have a
>> simple doubt.
>>
>> To extract I'm transforming my shape to the crs of the rasterstack to
>> get the same crs.
>>
>> What I'm concerned about is that before and after the changing in the
>> crs when I ask for the function proj4string I have a warning:
>>
>> "In proj4string(shape) : CRS object has comment, which is lost in output"
>>
>> I read that it's because of the development of the rgdal and that in
>> most cases there is no problem. But I would like confirmation that this
>> is not affecting my results of extraction. I believe that if the warning
>> came after the extraction, I should be more concerned. Could anyone
>> help?
>
> Not rgdal itself, but the underlying PROJ library, which has changed
> dramatically. The warnings should lead users to check at least three times
> before assuming that things are OK.
>
>>
>> The code Is like that:
>>
>> shape <- readOGR(dsn = ".", layer = "myshape")
>>
>> proj4string(shape)
>>
>> fn<-file.path("mypath\\cru_ts4.04.2011.2019.pre.dat.nc")
>>
>> rasbrick <- stack(fn)
>>
>
> I see for this publically available file, with current CRAN sp, rgdal and
> raster:
>
>> str(crs(rasbrick))
> Formal class 'CRS' [package "sp"] with 1 slot
>   ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
>
> but crucially:
>
>> cat(wkt(crs(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 WKT2-2019 representation also present in the CRS object suggests that
> the CRU_TS grid is read correctly.
>
>> shape2<- spTransform(shape, crs(rasbrick))
>>
>
> The success of this operation will depend on:
>
> cat(wkt(shape), "\n")
>
> being well-defined, the on there being an optimal path from the CRS of
> shape to that of rasbrick (up to PROJ 6, there was no choice in the
> operation, as the source and target CRS defined their relationship to
> WGS84). After spTransform(), you can use get_last_coordOp() to show the
> coordinate operation pipeline. In sp/rgdal, every effort has been made to
> prevent axis-swapping.
>
>> proj4string(shape2)
>
> Rather check with cat(wkt(shape2), "\n"), as proj4string() may issue a
> warning, as you are asking to look at the deprecated Proj4 representation.
>
>>
>> for (i in 1:length(rasbrick using layers)) {
>>  weath_dt[,2+i] = raster::extract(rasbrick[[i]], shape2, mean, na.rm=T)
>> }
>
> Because raster uses sp/rgdal, it may also issue warnings. The warnings are
> there to provoke questions like yours. Warnings may be muted for sp and
> rgdal by options("rgdal_show_exportToProj4_warnings"="none") before
> loading rgdal.
>
> Hope this clarifies,
>
> Roger
>
>>
>> Thanks.
>>
>> Pietro Andre Telatin Paschoalino
>> Doutorando em Ci�ncias Econ�micas da Universidade Estadual de Maring� - PCE.
>>
>>        [[alternative HTML version deleted]]
>>
>>
>
> --
> Roger Bivand
> 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
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