[R-sig-Geo] Simple doubt - CRS

[Roger Bivand]

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