Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Semantics of crop(SpatRaster, extent) on SpatRaster attributes #1302

Closed
rhgof opened this issue Oct 5, 2023 · 3 comments
Closed

Semantics of crop(SpatRaster, extent) on SpatRaster attributes #1302

rhgof opened this issue Oct 5, 2023 · 3 comments

Comments

@rhgof
Copy link

rhgof commented Oct 5, 2023

I am trying to work out the behavior of SpatRaster attributes when created from a netcdf file.

The example is below (the .nc file) that is a grid of readings off satellite
Before the crop I have all the data including the units.
After the crop - the attributes units, varname are dropped and can only be reconstructed using copy and assign functions.

Somewhat surprisingly when I create the raster in a different way, and crop, attribute units is preserved.
See the second example.

crop() clearly preserves some attributes, but not all and does not appear to do so consistently.

I am sure I am missing something here - as only a few weeks starting with Terra

I want to carry as many attributes forward as possible as part of developing a data pipeline to do multiple data source analysis on geo-graphic regions.

Many thanks

library(terra)
#> terra 1.7.46
theRast <- rast("/Volumes/Samples/InputData/cache/A.P1D.20230821T053000Z.aust.chl_gsm.nc")
theRast
#> class       : SpatRaster 
#> dimensions  : 7001, 10001, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 79.995, 180.005, -60.005, 10.005  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : A.P1D.20230821T053000Z.aust.chl_gsm.nc 
#> varname     : chl_gsm (Chlorophyll Concentration, GSM model) 
#> name        : chl_gsm 
#> unit        :  mg/m^3 
#> time (days) : 2023-08-21

ln <- longnames(theRast)
vn <- varnames(theRast)
un <- units(theRast)

theRast <- crop(theRast,ext(c(140,150,-50,-10)))
theRast
#> class       : SpatRaster 
#> dimensions  : 4000, 1000, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 140.005, 150.005, -49.995, -9.995  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> name        :      chl_gsm 
#> min value   :  0.001007809 
#> max value   : 98.747230530 
#> time (days) : 2023-08-21

# Tests
longnames(theRast) == ln
#> [1] FALSE
units(theRast) == un
#> [1] FALSE
varnames(theRast) == vn
#> [1] FALSE

longnames(theRast) <- ln
varnames(theRast) <- vn
units(theRast) <- un

# More Test
longnames(theRast) == ln
#> [1] TRUE
units(theRast) == un
#> [1] TRUE
varnames(theRast) == vn
#> [1] TRUE

theRast
#> class       : SpatRaster 
#> dimensions  : 4000, 1000, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 140.005, 150.005, -49.995, -9.995  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> varname     : chl_gsm (Chlorophyll Concentration, GSM model) 
#> name        :      chl_gsm 
#> min value   :  0.001007809 
#> max value   : 98.747230530 
#> unit        :       mg/m^3 
#> time (days) : 2023-08-21

Created on 2023-10-05 with reprex v2.0.2

R.version.string
#> [1] "R version 4.2.1 (2022-06-23)"
library(terra)
#> terra 1.7.46

r <- rast(nrows=180, ncols=360, nlyrs=1, xmin=-180, xmax=180, ymin=-90, ymax=90)
crs(r) <- "EPSG:7844"
r
#> class       : SpatRaster 
#> dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat GDA2020 (EPSG:7844)
set.names(r, "Test")
names(r)
#> [1] "Test"

units(r) <-"t/s"
units(r)
#> [1] "t/s"

varnames(r) <-"TestVar"
varnames(r)
#> [1] "TestVar"

longnames(r) <- "TestLong"
longnames(r)
#> [1] "TestLong"

rc <- crop(r,ext(c(140,150,-50,-10)))
rc
#> class       : SpatRaster 
#> dimensions  : 40, 10, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 140, 150, -50, -10  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat GDA2020 (EPSG:7844)
# EQUAL
names(rc) == names(r)
#> [1] TRUE
units(rc) == units(r)
#> [1] TRUE
crs(rc) == crs(r)
#> [1] TRUE
# NOT EQUAL
varnames(rc) == varnames(r)
#> [1] FALSE
longnames(rc) == longnames(r)
#> [1] FALSE

Created on 2023-10-05 with reprex v2.0.2

@rhijmans
Copy link
Member

rhijmans commented Oct 5, 2023

Thanks. The units getting lost was a bug and that is fixed now.
The var- and longnames are now also conserved if the input has a single data source (if there is only one var/longname).

@rhgof
Copy link
Author

rhgof commented Oct 6, 2023

Why are var/long-names not conserved if multiple layers and sources?
And why does the capture and re-assignment below not work?
Whereas the units can be re-assigned for multiple layers ok?

library(terra)
#> terra 1.7.46

theFiles = c("/Volumes/Samples/InputData/cache/A.P1D.20230821T053000Z.aust.chl_gsm.nc",
             "/Volumes/Samples/InputData/cache/A.P1D.20230822T053000Z.aust.chl_gsm.nc",
             "/Volumes/Samples/InputData/cache/A.P1D.20230823T053000Z.aust.chl_gsm.nc"
             )
theRast = rast(theFiles)
theRast
#> class       : SpatRaster 
#> dimensions  : 7001, 10001, 3  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 79.995, 180.005, -60.005, 10.005  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> sources     : A.P1D.20230821T053000Z.aust.chl_gsm.nc  
#>               A.P1D.20230822T053000Z.aust.chl_gsm.nc  
#>               A.P1D.20230823T053000Z.aust.chl_gsm.nc  
#> varnames    : chl_gsm (Chlorophyll Concentration, GSM model) 
#>               chl_gsm (Chlorophyll Concentration, GSM model) 
#>               chl_gsm (Chlorophyll Concentration, GSM model) 
#> names       : chl_gsm, chl_gsm, chl_gsm 
#> unit        :  mg/m^3,  mg/m^3,  mg/m^3 
#> time (days) : 2023-08-21 to 2023-08-23
theUnits = units(theRast)
longNames <- longnames(theRast)
vn <-varnames(theRast)

theRast<-crop(theRast,ext(c(140,150,-50,-10)))
theRast
#> class       : SpatRaster 
#> dimensions  : 4000, 1000, 3  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 140.005, 150.005, -49.995, -9.995  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> names       :      chl_gsm,      chl_gsm,      chl_gsm 
#> min values  :  0.001007809,  0.001007874,  0.001009486 
#> max values  : 98.747230530, 95.338264465, 98.991081238 
#> time (days) : 2023-08-21 to 2023-08-23

varnames(theRast) <- vn
#> Error: [varnames<-,SpatRaster] cannot set these names
longnames(theRast) <- longNames
#> Error: [longnames<-,SpatRaster] cannot set these names
units(theRast) <- theUnits

theRast
#> class       : SpatRaster 
#> dimensions  : 4000, 1000, 3  (nrow, ncol, nlyr)
#> resolution  : 0.01, 0.01  (x, y)
#> extent      : 140.005, 150.005, -49.995, -9.995  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> names       :      chl_gsm,      chl_gsm,      chl_gsm 
#> min values  :  0.001007809,  0.001007874,  0.001009486 
#> max values  : 98.747230530, 95.338264465, 98.991081238 
#> unit        :       mg/m^3,       mg/m^3,       mg/m^3 
#> time (days) : 2023-08-21 to 2023-08-23

Created on 2023-10-06 with reprex v2.0.2

rhijmans added a commit that referenced this issue Oct 6, 2023
rhijmans added a commit that referenced this issue Oct 6, 2023
@rhijmans
Copy link
Member

rhijmans commented Oct 6, 2023

var/long-names are properties of "data sources". These may exist in the input, but since the output is always a single data source, they cannot always be maintained in the output. "terra" now keeps them if they are all the same, as in your example. In other cases you would need to process variable by variable (data source); perhaps using SpatRasterDatasets. I have also expanded the docs to clarify this a bit.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants