R Under development (unstable) (2023-08-09 r84924 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > Sys.setenv(TZ="UTC") > suppressPackageStartupMessages(library(stars)) > set.seed(13521) # runif > tif = system.file("tif/L7_ETMs.tif", package = "stars") > (x_ = read_stars(c(tif,tif))) # FIXME: not what you'd expect stars object with 3 dimensions and 2 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 L7_ETMs.tif.1 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > (x = read_stars(tif)) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > # image(x) > gdal_crs(tif) Coordinate Reference System: User input: PROJCS["SIRGAS 2000 / UTM zone 25S",GEOGCS["SIRGAS 2000",DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6674"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4674"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","31985"]] wkt: PROJCRS["SIRGAS 2000 / UTM zone 25S", BASEGEOGCRS["SIRGAS 2000", DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4674]], CONVERSION["UTM zone 25S", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-33, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9996, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",500000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",10000000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1]], ID["EPSG",31985]] > plot(x) > plot(x, join_zlim = FALSE) > x %>% st_set_dimensions(names = c('a', 'b', 'c')) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x/y a 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] b 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] c 1 6 NA NA NA NA > st_get_dimension_values(x, 3) [1] 1 2 3 4 5 6 > > (x1 = st_set_dimensions(x, "band", values = c(1,2,3,4,5,7), names = "band_number", point = TRUE)) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point values x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [y] band_number 1 6 NA NA NA TRUE 1,...,7 > rbind(c(0.45,0.515), c(0.525,0.605), c(0.63,0.69), c(0.775,0.90), c(1.55,1.75), c(2.08,2.35)) %>% + units::set_units(um) -> bw # units::set_units(µm) -> bw > # set bandwidth midpoint: > (x2 = st_set_dimensions(x, "band", values = 0.5 * (bw[,1]+bw[,2]), + names = "bandwidth_midpoint", point = TRUE)) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE bandwidth_midpoint 1 6 NA NA udunits TRUE values x/y x NULL [x] y NULL [y] bandwidth_midpoint 0.4825 [um],...,2.215 [um] > # set bandwidth intervals: > (x3 = st_set_dimensions(x, "band", values = make_intervals(bw), names = "bandwidth")) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE bandwidth 1 6 NA NA udunits NA values x/y x NULL [x] y NULL [y] bandwidth [0.45,0.515) [um],...,[2.08,2.35) [um] > > x + x stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 2 108 138 137.8248 172 510 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > x * x stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 2916 4761 5512.41 7396 65025 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > x[,,,1:3] stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 21 58 70 70.36041 83 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 3 NA NA NA NA > x[,1:100,100:200,] stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 13 54 65 67.21531 77 252 dimension(s): from to offset delta refsys point x/y x 1 100 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 100 200 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > sqrt(x) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 7.348469 8.306624 8.094006 9.273618 15.96872 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > st_apply(x, 3, min) stars object with 1 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. min 1 3 15 18.5 29.25 47 dimension(s): from to band 1 6 > st_apply(x, 1:2, max) stars object with 2 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. max 55 85 96 99.36018 113 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] > st_apply(x, 1:2, range) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 50 71 72.43565 96 255 dimension(s): from to offset delta refsys point x/y range 1 2 NA NA NA NA x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] > > geomatrix = system.file("tif/geomatrix.tif", package = "stars") > x = read_stars(geomatrix) > y = st_transform(x, st_crs(4326)) > st_coordinates(x)[1:10,] x y 1 1841000 1144000 2 1841002 1143995 3 1841003 1143990 4 1841005 1143985 5 1841006 1143980 6 1841008 1143975 7 1841009 1143970 8 1841011 1143965 9 1841012 1143960 10 1841014 1143955 > > nc = system.file("nc/tos_O1_2001-2002.nc", package = "stars") > if (nc != "" && require(PCICt, quietly = TRUE)) { + print(x <- read_stars(nc)) + print(st_bbox(x)) + s = st_as_stars(st_bbox(x)) # inside = NA + print(st_bbox(s)) + s = st_as_stars(st_bbox(x), inside = TRUE) + print(st_bbox(s)) + s = st_as_stars(st_bbox(x), inside = FALSE) + print(st_bbox(s)) + (s = st_as_stars(st_bbox(x), dy = 1)) + print(st_bbox(s)) + print(identical(st_as_stars(st_bbox(x), dx = 1), st_as_stars(st_bbox(x), dy = 1))) + s = st_as_stars(st_bbox(x), dx = 10) + print(st_bbox(s)) + s = st_as_stars(st_bbox(x), dx = 20) + print(st_bbox(s)) + x1 = x + st_crs(x1) = "OGC:CRS84" + print(identical(st_as_stars(st_bbox(x1), dx = 1), st_as_stars(st_bbox(x1), dx = units::set_units(1, degree)))) + + df = as.data.frame(x) + if (require(units, quietly = TRUE)) + print(units::drop_units(x)) + + print(dimnames(x)) + dimnames(x) <- letters[1:3] + print(dimnames(x)) + } # PCICt > > st_as_stars() stars object with 2 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. values 0 0 0 0 0 0 dimension(s): from to offset delta refsys x/y x 1 360 -180 1 WGS 84 [x] y 1 180 90 -1 WGS 84 [y] > > # multiple sub-datasets: > nc_red = system.file("nc/reduced.nc", package = "stars") > red = read_stars(nc_red) sst, anom, err, ice, > red stars object with 4 dimensions and 4 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's sst [°C] -1.80 -0.03 13.655 12.9940841 24.8125 32.97 4448 anom [°C] -10.16 -0.58 -0.080 -0.1855812 0.2100 2.99 4448 err [°C] 0.11 0.16 0.270 0.2626872 0.3200 0.84 4448 ice [percent] 0.01 0.47 0.920 0.7178118 0.9600 1.00 13266 dimension(s): from to offset delta refsys x/y x 1 180 -1 2 NA [x] y 1 90 90 -2 NA [y] zlev 1 1 0 [m] NA NA time 1 1 1981-12-31 UTC NA POSIXct > plot(red) > > x = st_xy2sfc(read_stars(tif)[,1:10,1:10,], as_points = FALSE) > st_bbox(x) xmin ymin xmax ymax 288776.3 9120475.8 289061.3 9120760.8 > x = read_stars(tif) > merge(split(x, "band")) stars object with 3 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. X1.X2.X3.X4.X5.X6 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point values x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL attributes 1 6 NA NA NA NA X1,...,X6 x/y x [x] y [y] attributes > > read_stars(c(tif,tif)) # merges as attributes stars object with 3 dimensions and 2 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 1 54 69 68.91242 86 255 L7_ETMs.tif.1 1 54 69 68.91242 86 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > read_stars(c(tif,tif), along = "sensor") stars object with 4 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 47 65 76 77.3419 87 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA sensor 1 2 NA NA NA NA > read_stars(c(tif,tif), along = 4) stars object with 4 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 47 65 76 77.3419 87 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA new_dim 1 2 NA NA NA NA > read_stars(c(tif,tif), along = "band") stars object with 3 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 47 65 76 77.3419 87 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 12 NA NA NA NA > read_stars(c(tif,tif), along = 3) stars object with 3 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: Min. 1st Qu. Median Mean 3rd Qu. Max. L7_ETMs.tif 47 65 76 77.3419 87 255 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 12 NA NA NA NA > > # cut: > tif = system.file("tif/L7_ETMs.tif", package = "stars") > x = read_stars(tif) > cut(x, c(0, 50, 100, 255)) stars object with 3 dimensions and 1 attribute attribute(s): L7_ETMs.tif (0,50] :156060 (50,100] :503764 (100,255]: 77264 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] band 1 6 NA NA NA NA > cut(x[,,,1,drop=TRUE], c(0, 50, 100, 255)) stars object with 2 dimensions and 1 attribute attribute(s): L7_ETMs.tif (0,50] : 1 (50,100] :117134 (100,255]: 5713 dimension(s): from to offset delta refsys point x/y x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x] y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y] > plot(cut(x[,,,1,drop=TRUE], c(0, 50, 100, 255))) > > st_bbox(st_dimensions(x)) xmin ymin xmax ymax 288776.3 9110728.8 298722.8 9120760.8 > x[x < 0] = NA > x[is.na(x)] = 0 > > # c: > f = system.file("netcdf/avhrr-only-v2.19810902.nc", package = "starsdata") > if (FALSE && f != "") { + files = c("avhrr-only-v2.19810901.nc", + "avhrr-only-v2.19810902.nc", + "avhrr-only-v2.19810903.nc", + "avhrr-only-v2.19810904.nc", + "avhrr-only-v2.19810905.nc", + "avhrr-only-v2.19810906.nc", + "avhrr-only-v2.19810907.nc", + "avhrr-only-v2.19810908.nc", + "avhrr-only-v2.19810909.nc") + l = list() + for (f in files) { + from = system.file(paste0("netcdf/", f), package = "starsdata") + l[[f]] = read_stars(from, sub = c("sst", "anom")) + } + ret = do.call(c, l) + print(ret) + ret = adrop(c(l[[1]], l[[2]], l[[3]], along = list(times = as.Date("1981-09-01") + 0:2))) + print(ret) + #ret = adrop(adrop(c(l[[1]], l[[2]], l[[3]], along = "times"))) + #print(ret) + } > > st_dimensions(list(matrix(1, 4, 4))) # st_dimensions.default from to point values X1 1 1 FALSE 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 > > if (FALSE && require("starsdata", quietly = TRUE)) { + # curvilinear: + s5p = system.file( + "sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", + package = "starsdata") + print(s5p) + lat_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/latitude") + lon_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/longitude") + nit_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column") + lat = read_stars(lat_ds) + lon = read_stars(lon_ds) + nit = read_stars(nit_ds) + nit[[1]][nit[[1]] > 9e+36] = NA + + ll = setNames(c(lon, lat), c("x", "y")) + nit.c = st_as_stars(nit, curvilinear = ll) + print(nit.c) + + s5p = system.file( + "sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", + package = "starsdata") + nit.c2 = read_stars(s5p, + sub = "//PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column", + curvilinear = c("//PRODUCT/latitude", "//PRODUCT/longitude")) + print(all.equal(nit.c, nit.c2)) + } > > proc.time() user system elapsed 5.25 0.59 5.84