Simple doubt - CRS

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Simple doubt - CRS

Pietro Andre Telatin Paschoalino
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Simple doubt - CRS

Roger Bivand
Administrator
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@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: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Simple doubt - CRS

Pietro Andre Telatin Paschoalino
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Simple doubt - CRS

Roger Bivand
Administrator
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 <[hidden email]>
> Enviado: quarta-feira, 24 de março de 2021 11:11
> Para: Pietro Andre Telatin Paschoalino <[hidden email]>
> Cc: [hidden email] <[hidden email]>
> 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@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: [hidden email]
> 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: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway