Simple doubt - CRS

Enviado: quarta-feira, 24 de mar?o de 2021 11:11
Para: Pietro Andre Telatin Paschoalino <Pietro_telato at hotmail.com>
Cc: r-sig-geo at r-project.org <r-sig-geo at 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 at 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 at nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en