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.

netbsd-srcmastr pushed a commit to NetBSD/pkgsrc that referenced this issue Dec 10, 2024
# version 1.7-83

## bug fixes

- `flip(direction="vertical")` failed in some cases
  [#1518](rspatial/terra#1518) by Ed Carnell

- `zonal(as.raster=TRUE)` failed when the zonal raster was categorical
  [1514](rspatial/terra#1514) by Jessi L
  Brown

- `distance<data.frame,data.frame>` and `<matrix,matrix>` ignored the
  unit
  argument. [#1545](rspatial/terra#1545) by
  Wencheng Lau-Medrano

- NetCDF files with month time-step encode from 0-11 made R crash
  [#1544](rspatial/terra#1544) by Martin
  Holdrege

- `split<SpatVector>` only worked well if the split field was of type
  character. [#1530](rspatial/terra#1530) by
  Igor Graczykowski

- `gridDist` (and probably some other methods) emitted a "cannot
  overwrite existing file" error when processing large datasets
  [#1522](rspatial/terra#1522) by Clare
  Pearson

- `terrain` did not accept multiple variables
  [#1561](rspatial/terra#1561) by Michael
  Mahoney

- `rotate` was vulnerable to an integer overflow
  [#1562](rspatial/terra#1562) by Sacha
  Ruzzante

- `getTileExtents` could return overlapping tiles or tiles with gaps
  due to floating point
  imprecision. [#1564](rspatial/terra#1564)
  by Michael Sumner


## enhancements

- `as.list<SpatRasterDataset>` sets the names of the list
  [#1513](rspatial/terra#1513)

- a SpatVectorCollection can now be subset with its names; and if made
  from a list it takes the names from the list.
  [1515](rspatial/terra#1515) by jedgroev

- argument `fill_range` to plot<SpatRaster> and `plot<SpatVector>` to
  use the color of the extreme values of the specified range
  [#1553](rspatial/terra#1553) by Mike
  Koontz

- plet<SpatRaster> can now handle rasters with a "local" (Cartesian)
  CRS. [#1570](rspatial/terra#1570) by
  Augustin Lobo.

## new

- `map-region` returns the coordinates of the axes position of a map
  created with `plot<Spat*>`
  [https://github.com/rspatial/terra/issues/1517](https://github.com/rspatial/terra/issues/1517)
  by Daniel Schuch

- `polys<leaflet>` method
  [#1543](rspatial/terra#1543) by Márcia
  Barbosa

- `plot<SpatVectorCollection>` method
  [#1532](rspatial/terra#1532) by jedgroev

- `add_mtext` to add text around the margins of a
  map. [#1567](rspatial/terra#1567) by
  Daniel Schuch

# version 1.7-78

Released 2023-05-22

## bug fixes

- `writeVector` and `readVector` better handle empty geopackage layers
  [#1426](rspatial/terra#1426) by Andrew
  Gene Brown.

- `writeCDF` only wrote global variables if there was more than one
  [#1443](rspatial/terra#1443) by Daniel
  Schlaepfer

- `rasterize` with "by" returned odd layernames
  [#1435](rspatial/terra#1435) by Philippe
  Massicotte

- `convHull`, `minCircle` and `minRect` with a zero-row SpatVector
  crashed R [#1445](rspatial/terra#1445) by
  Andrew Gene Brown

- `rangeFill` with argument `circular=TRUE` did not work properly
  [#1460](rspatial/terra#1460) by Alice

- `crs(describe = TRUE)` returned an mis-ordered extent
  [#1485](rspatial/terra#1485) by Dimitri
  Falk

- `tapp` with a custom function and an index like "yearmonths" could
  shift time for not considering the time
  zone. [#1483](rspatial/terra#1483) by Finn
  Roberts

- `plot<SpatRaster>` could fail when there were multiple values with
  very small differences
  [#1491](rspatial/terra#1491) by srfall

- `as.data.frame<SpatRaster>` with "xy=TRUE" and "wide=FALSE" could
  fail if coordinates were very similar
  [#1476](rspatial/terra#1476) by Pascal
  Oettli

- `rasterizeGeom` now returns the correct layer name
  [#1472](rspatial/terra#1472) by
  HRodenhizer

- `cellSize` with "mask=TRUE" failed if the output was to be written
  to a temp file
  [#1496](rspatial/terra#1496) by Pascal
  Sauer

- `ext<SpatVectorProxy>` did not return the full extent
  [#1501](rspatial/terra#1501) by
  erkent-carb


## enhancements

- `extract` has new argument "small=TRUE" to allow for strict use of
  "touches=FALSE"
  [#1419](rspatial/terra#1419) by Floris
  Vanderhaeghe.

- `as.list<SpatRaster>` has new argument "geom=NULL"

- `rast<list>` now recognizes (x, y, z) base R "image" structures
  [stackoverflow]
  (https://stackoverflow.com/questions/77949551/rspatial-convert-a-grid-list-to-a-raster-using-terra)
  by Ignacio Marzan.

- `inset` has new arguments "offset" and "add"
  [#1422](rspatial/terra#1422) by Armand-CT

- `expanse<SpatRaster>` has argument `usenames`
  [#1446](rspatial/terra#1446) by Bappa Das

- the default color palette is now `terra::map.pal("viridis")` instead
  of `terrain.colors`. The default can be changes with
  `options(terra.pal=...)`
  [#1474](rspatial/terra#1474) by Derek
  Friend

- `as.list<SpatRasterDataset>` now returns a named
  list. [#1513](rspatial/terra#1513) by Eric
  R. Scott


## new

- `bestMatch<SpatRaster>` method

- argument "pairs=TRUE" to `cells` [https://github.com/rspatial/terra/issues/1487](https://github.com/rspatial/terra/issues/1487) by Floris Vanderhaeghe

- `add_grid` to add a grid to a map


# version 1.7-71

Released 2023-01-31

## bug fixes

- k_means did not work if there were NAs
  [#1314](rspatial/terra#1314) by Jakub
  Nowosad

- `layerCor` with a custom function did not work anymore
  [#1387](rspatial/terra#1387) by Jakub
  Nowosad

- `plet` broke when using "panel=TRUE"
  [#1384](rspatial/terra#1384) by Elise
  Hellwig

- using /vis3/ to open a SpatRaster did not work
  [#1382](rspatial/terra#1382) by Mike
  Koontz

- `plot<SpatRaster>(add=TRUE)` sampled the raster data without
  considering the extent of the
  map. [#1394](rspatial/terra#1394) by
  Márcia Barbosa

- `plot<SpatRaster>(add=TRUE)` now only considers the first layer of a
  multi-layer SpatRaster
  [1395](rspatial/terra#1395) by Márcia
  Barbosa

- `set.cats` failed with a tibble was used instead of a data.frame
  [#1406](rspatial/terra#1406) by Mike
  Koontz

- `polys` argument "alpha" was ignored if a single color was
  used. [#1413](rspatial/terra#1413) by
  Derek Friend

- `query` ignore the "vars" argument if all rows were
  selected. [#1398](rspatial/terra#1398) by
  erkent-carb.

- `spatSample` ignored "replace=TRUE" with random sampling,
  na.rm=TRUE, and a sample size larger than the non NA
  cells. [#1411](rspatial/terra#1411) by
  Babak Naimi

- `spatSample` sometimes returned fewer values than requested and
  available for lonlat
  rasters. [#1396](rspatial/terra#1396) by
  Márcia Barbosa.


## enhancements

- `vect<character>` now has argument "opts" for GDAL open options,
  e.g. to declare a file
  encoding. [#1389](rspatial/terra#1389) by
  Mats Blomqvist

- `plot(plg=list(tic=""))` now allows choosing alternative continuous
  legend tic-mark styles ("in", "out", "through" or "none")

- `makeTiles` has new argument "buffer"
  [#1408](rspatial/terra#1408) by Joy
  Flowers.


## new

- `prcomp<SpatRaster>` method
  [#1361](rspatial/terra#1361 (comment))
  by Jakub Nowosad

- `add_box` to add a box around the map. The box is drawn where the
  axes are, not around the plotting region.

- `getTileExtents` provides the extents of tiles. These may be used in
parallelization. See [#1391](https://github.com/rspa
tial/terra/issues/1391) by Alex Ilich.


# version 1.7-65

Released 2023-12-15

## bug fixes

- `flip` with argument `direction="vertical"` filed in some cases with
   large rasters processed in chunks
   [0b714b0](rspatial/terra@0b714b0)
   by Dulci on [stackoveflow](
   https://stackoverflow.com/questions/77304534/rspatial-terraflip-error-when-flipping-a-multi-layer-spatrast-object)

- SpatRaster now correctly handles `NA & FALSE` and `NA | TRUE`
  [#1316](rspatial/terra#1316) by John Baums

- `set.names` wasn't working properly for SpatRasterDataset or
  SpatRasterCollection
  [#1333](rspatial/terra#1333) by Derek Friend

- `extract` with argument "layer" not NULL shifted the layers
  [#1332](rspatial/terra#1332) by Ewan
  Wakefield

- `terraOptions` did not capture "memmin" on
  -[stackoverflow](https://stackoverflow.com/questions/77552234/controlling-chunk
  -size-in-terra) by dww

- `rasterize` with points and a built-in function could crash if no
  field was used
  [#1369](rspatial/terra#1369) by
  anjelinejeline


## enhancements

- `mosaic` can now use `fun="modal"`

- `rast<matrix> and rast<data.frame>` now have option 'type="xylz"
  [#1318](rspatial/terra#1318) by Agustin
  Lobo

- `extract<SpatRaster,SpatVector>` can now use multiple summarizing
  functions [#1335](rspatial/terra#1335) by
  Derek Friend

- `disagg` and `focal` have more optimistic memory requirement
  estimation [#1334](rspatial/terra#1334) by
  Mikko Kuronen

## new

- `k_means<SpatRaster>` method
  [#1314](rspatial/terra#1314) by Agustin
  Lobo

- `princomp<SpatRaster>` method
  [#1361](rspatial/terra#1361) by Alex Ilich

- `has.time<SpatRaster>` method

- new argument "raw=FALSE" to `rast`, `sds`, and `sprc` to allow
  ignoring scale and offset
  [1354](rspatial/terra#1354) by Insang Song


# version 1.7-55

Released 2023-10-14

## bug fixes

- `mosaic` ignored the filename argument if the SpatRasterCollection
  only had a single SpatRaster
  [#1267](rspatial/terra#1267) by Michael
  Mahoney

- Attempting to use `extract` with a raster file that had been deleted
  crashed R. [#1268](rspatial/terra#1268) by
  Derek Friend

- `split<SpatVector,SpatVector>` did not work well in all
  cases. [#1256](rspatial/terra#1256) by
  Derek Corcoran Barrios

- `intersect` with two SpatVectors crashed R if there was a date/time
variable [#1273]( rspatial/terra#1273) by
Dave Dixon

- "values=FALSE" was ignored by
  `spatSample<SpatRaster>(method="weights")`
  [#1275](rspatial/terra#1275) by François
  Rousseu

- `coltab<-` again works with a list as value
[#1280](rspatial/terra#1280) by Diego
Hernangómez

- `stretch` with histogram equalization was not memory-safe
  [#1305](rspatial/terra#1305) by Evan Hersh

- `plot` now resets the "mar" parameter
  [#1297](rspatial/terra#1297) by Márcia
  Barbosa

- `plotRGB` ignored the "smooth" argument
  [#1307](rspatial/terra#1307) by Timothée
  Giraud


## enhancements

- argument "gdal" in `project` was renamed to "use_gdal"
  [#1269](rspatial/terra#1269) by Stuart
  Brown.

- SpatVector attributes can now be stored as an ordered factor
  [#1277](rspatial/terra#1277) by Ben Notkin

- `plot<SpatVector>` now uses an "interval" legend when breaks are
  supplied [#1303](rspatial/terra#1303) by
  Gonzalo Rizzo

- `crop<SpatRaster>` now keeps more metadata, including variable names
  [#1302](rspatial/terra#1302) by rhgof

- `extract(fun="table")` now returns an easier to use data.frame
[#1294](rspatial/terra#1294) by Fernando
Aramburu.


## new
- `metags<-` and `metags` to set arbitrary SpatRaster/file level
   metadata [#1304](https://github.com/rspatial/terra/issues/ 1304) by
   Francesco Chianucci

# version 1.7-46

Released 2023-09-06

## bug fixes

- `plot<SpatVector>` used the wrong main label in some cases
  [#1210](rspatial/terra#1210) by Márcia
  Barbosa

- `plotRGB` failed with an "ext=" argument
  [#1228](rspatial/terra#1228) by Dave Edge

- `rast<array>` failed badly when the array had less than three
  dimensions. [#1254](rspatial/terra#1254)
  by andreimirt.

- `all.equal` for a SpatRaster with multiple layers
[#1236](rspatial/terra#1236) by Sarah
Endicot t

- `zonal(wide=FALSE)` could give wrong results if the zonal SpatRaster
  had "layer" as
  layername. [#1251](rspatial/terra#1251) by
  Jeff Hanson

- `panel` now support argument "range"
  [#141](rspatial/terra#1241) by Jakub
  Nowosad

- `rasterize` with `by=` returned wrong layernames if the by field was
  not sorted [#1266](rspatial/terra#1266) by
  Sebastian Dunnett

- `mosaic` with multiple layers was not correct
  [#1262](rspatial/terra#1262) by
  Jean-Romain


## enhancements

- `wrap<SpatRaster>` now stores color tables
  [#1215](rspatial/terra#1215) by Patrick
  Brown

- `global` now has a "maxcell" argument
  [#1213](rspatial/terra#1213) by Alex Ilich

- `layerCor` with fun='pearson' now returns output with the layer
  names [#1206](rspatial/terra#1206)

- `vrt` now has argument "set_names"
  [#1244](rspatial/terra#1244) by sam-a-levy

- `vrt` now has argument "return_filename"
  [#1258](rspatial/terra#1258) by Krzysztof
  Dyba

- `project<SpatRaster>` has new argument "by_util" exposing the GDAL
  warp utility [#1222](rspatial/terra#1222) by
  Michael Sumner.


## new
- `compareGeom` for list and SpatRasterCollection
  [#1207](rspatial/terra#1207) by Sarah
  Endicott

- `is.rotated<SpatRaster>` method
  [#1229](rspatial/terra#1229) by Andy Lyons

- `forceCCW<SpatVector>` method to force counter-clockwise orientation
  of polygons [#1249](rspatial/terra#1249)
  by srfall.

- `vrt_tiles` returns the filenames of the tiles in a vrt file
  [#1261](rspatial/terra#1261) by Derek
  Friend

- `extractAlong` to extract raster cell values for a line that are
  ordered along the
  line. [#1257](rspatial/terra#1257) by
  adamkc.
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