R Under development (unstable) (2023-07-09 r84667 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. > suppressPackageStartupMessages(library(sf)) > > tif = system.file("tif/geomatrix.tif", package = "sf") > > gdal_metadata(tif) [1] "AREA_OR_POINT=Point" > gdal_metadata(tif, NA_character_) [1] "IMAGE_STRUCTURE" "DERIVED_SUBDATASETS" "" > try(gdal_metadata(tif, "wrongDomain")) Error in gdal_metadata(tif, "wrongDomain") : domain_item[1] not found in available metadata domains > gdal_metadata(tif, c("IMAGE_STRUCTURE")) $INTERLEAVE [1] "BAND" attr(,"class") [1] "gdal_metadata" > try(length(gdal_metadata(tif, c("DERIVED_SUBDATASETS")))) # fails on Fedora 26 [1] 2 > > if (require(stars, quietly = TRUE)) { + tif = system.file("tif/geomatrix.tif", package = "sf") + r = read_stars(tif) + d = (st_dimensions(r)) + gt = c(1841001.75, 1.5, -5, 1144003.25, -5, -1.5) + x1 = st_as_sfc(d, as_points = TRUE, use_cpp = TRUE, geotransform = gt) + x2 = st_as_sfc(d, as_points = TRUE, use_cpp = FALSE, geotransform = gt) + print(identical(x1, x2)) + y1 = st_as_sfc(d, as_points = FALSE, use_cpp = TRUE, geotransform = gt) + y2 = st_as_sfc(d, as_points = FALSE, use_cpp = FALSE, geotransform = gt) + print(identical(y1, y2)) + + # rectilinear grid: + m = matrix(1:20, nrow = 5, ncol = 4) + x = c(0,0.5,1,2,4,5) + y = c(0.3,0.5,1,2,2.2) + r = st_as_stars(list(m = m), dimensions = st_dimensions(x = x, y = y, .raster = c("x", "y"))) + print(st_as_sfc(st_dimensions(r), as_points = TRUE)) + print(st_as_sfc(st_dimensions(r), as_points = FALSE)) + + # curvilinear grid: + lon = st_as_stars(matrix(1:5, 4, 5, byrow = TRUE)) + lat = st_as_stars(matrix(1:4, 4, 5)) + ll = c(X1 = lon, X2 = lat) + curv = st_as_stars(st_as_stars(t(m)), curvilinear = setNames(ll, c("X1", "X2"))) + print(st_as_sfc(st_dimensions(curv), as_points = TRUE)) + print(st_as_sfc(st_dimensions(curv), as_points = FALSE)) + + demo(nc, echo = FALSE, ask = FALSE) + print(x <- st_rasterize(nc)) # default grid: + print(p <- st_as_sf(x, as_points = FALSE)) # polygonize: follow raster boundaries + print(p <- st_as_sf(x, as_points = FALSE, use_integer = TRUE)) # polygonize integers: follow raster boundaries + print(try(p <- st_as_sf(x, as_points = TRUE))) # polygonize: contour, requies GDAL >= 2.4.0 + if (utils::packageVersion("stars") >= "0.2-1") { + write_stars(read_stars(tif), tempfile(fileext = ".tif")) + write_stars(read_stars(tif, proxy = TRUE), tempfile(fileext = ".tif")) + write_stars(read_stars(tif, proxy = TRUE), tempfile(fileext = ".tif"), chunk_size = c(200,200)) + na.tif = read_stars(system.file("tif/na.tif", package = "stars")) + write_stars(na.tif, "na.tif") + write_stars(na.tif, "na.tif", NA_value = -999) + na.tif = read_stars(system.file("tif/na.tif", package = "stars"), NA_value = -999) + write_stars(na.tif, "na.tif") + write_stars(na.tif, "na.tif", NA_value = -999) + na.tif = read_stars(system.file("tif/na.tif", package = "stars"), NA_value = -999, proxy = TRUE) + write_stars(na.tif, "na.tif") + write_stars(na.tif, "na.tif", NA_value = -999) + } + # https://github.com/mtennekes/tmap/issues/368 + if (utils::packageVersion("stars") > "0.4-0") { + lc = system.file('tif/lc.tif', package = 'stars') + if (lc != "") { + r = read_stars(lc, RAT = "Land Cover Class") + r <- droplevels(r) + } + } + } [1] TRUE [1] TRUE Geometry set for 20 features Geometry type: POINT Dimension: XY Bounding box: xmin: 0.25 ymin: 0.4 xmax: 4.5 ymax: 2.1 CRS: NA First 5 geometries: POINT (0.25 0.4) POINT (0.75 0.4) POINT (1.5 0.4) POINT (3 0.4) POINT (4.5 0.4) Geometry set for 20 features Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0 ymin: 0.3 xmax: 5 ymax: 2.2 CRS: NA First 5 geometries: POLYGON ((0 0.3, 0.5 0.3, 0.5 0.5, 0 0.5, 0 0.3)) POLYGON ((0.5 0.3, 1 0.3, 1 0.5, 0.5 0.5, 0.5 0... POLYGON ((1 0.3, 2 0.3, 2 0.5, 1 0.5, 1 0.3)) POLYGON ((2 0.3, 4 0.3, 4 0.5, 2 0.5, 2 0.3)) POLYGON ((4 0.3, 5 0.3, 5 0.5, 4 0.5, 4 0.3)) Geometry set for 20 features Geometry type: POINT Dimension: XY Bounding box: xmin: 1 ymin: 1 xmax: 5 ymax: 4 Geodetic CRS: WGS 84 First 5 geometries: POINT (1 1) POINT (1 2) POINT (1 3) POINT (1 4) POINT (2 1) Geometry set for 20 features Geometry type: POLYGON Dimension: XY Bounding box: xmin: 0.5 ymin: 0.5 xmax: 5.5 ymax: 4.5 Geodetic CRS: WGS 84 First 5 geometries: POLYGON ((0.5 0.5, 0.5 1.5, 1.5 1.5, 1.5 0.5, 0... POLYGON ((0.5 1.5, 0.5 2.5, 1.5 2.5, 1.5 1.5, 0... POLYGON ((0.5 2.5, 0.5 3.5, 1.5 3.5, 1.5 2.5, 0... POLYGON ((0.5 3.5, 0.5 4.5, 1.5 4.5, 1.5 3.5, 0... POLYGON ((1.5 0.5, 1.5 1.5, 2.5 1.5, 2.5 0.5, 1... stars object with 2 dimensions and 12 attributes attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's AREA 0.042 0.108 0.142 1.451932e-01 0.181 0.241 30904 PERIMETER 0.999 1.461 1.716 1.786110e+00 2.004 3.640 30904 CNTY_ 1825.000 1907.000 1989.000 1.998403e+03 2085.000 2241.000 30904 CNTY_ID 1825.000 1907.000 1989.000 1.998403e+03 2085.000 2241.000 30904 FIPSNO 37001.000 37049.000 37101.000 3.710042e+04 37153.000 37199.000 30904 CRESS_ID 1.000 25.000 51.000 5.071206e+01 77.000 100.000 30904 BIR74 248.000 1323.000 2648.000 3.791637e+03 4139.000 21588.000 30904 SID74 0.000 3.000 5.000 7.891985e+00 10.000 44.000 30904 NWBIR74 1.000 297.000 844.000 1.246210e+03 1396.000 8027.000 30904 BIR79 319.000 1606.000 3108.000 4.852046e+03 5400.000 30757.000 30904 SID79 0.000 3.000 6.000 9.584098e+00 13.000 57.000 30904 NWBIR79 3.000 360.000 1058.000 1.604642e+03 1524.000 11631.000 30904 dimension(s): from to offset delta refsys point x/y x 1 461 -84.3239 0.0192484 NAD27 FALSE [x] y 1 141 36.5896 -0.0192484 NAD27 FALSE [y] Simple feature collection with 34097 features and 12 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: -84.32385 ymin: 33.87563 xmax: -75.45034 ymax: 36.58965 Geodetic CRS: NAD27 First 10 features: AREA PERIMETER CNTY_ CNTY_ID FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 1 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 2 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 3 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 4 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 5 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 6 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 7 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 8 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 9 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 10 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 SID79 NWBIR79 geometry 1 0 19 POLYGON ((-81.66757 36.5896... 2 0 19 POLYGON ((-81.64833 36.5896... 3 0 19 POLYGON ((-81.62908 36.5896... 4 0 19 POLYGON ((-81.60983 36.5896... 5 0 19 POLYGON ((-81.59058 36.5896... 6 0 19 POLYGON ((-81.57133 36.5896... 7 0 19 POLYGON ((-81.55208 36.5896... 8 0 19 POLYGON ((-81.53283 36.5896... 9 0 19 POLYGON ((-81.51359 36.5896... 10 0 19 POLYGON ((-81.49434 36.5896... Simple feature collection with 34097 features and 12 fields Geometry type: POLYGON Dimension: XY Bounding box: xmin: -84.32385 ymin: 33.87563 xmax: -75.45034 ymax: 36.58965 Geodetic CRS: NAD27 First 10 features: AREA PERIMETER CNTY_ CNTY_ID FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 1 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 2 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 3 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 4 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 5 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 6 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 7 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 8 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 9 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 10 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 SID79 NWBIR79 geometry 1 0 19 POLYGON ((-81.66757 36.5896... 2 0 19 POLYGON ((-81.64833 36.5896... 3 0 19 POLYGON ((-81.62908 36.5896... 4 0 19 POLYGON ((-81.60983 36.5896... 5 0 19 POLYGON ((-81.59058 36.5896... 6 0 19 POLYGON ((-81.57133 36.5896... 7 0 19 POLYGON ((-81.55208 36.5896... 8 0 19 POLYGON ((-81.53283 36.5896... 9 0 19 POLYGON ((-81.51359 36.5896... 10 0 19 POLYGON ((-81.49434 36.5896... Simple feature collection with 34097 features and 12 fields Geometry type: POINT Dimension: XY Bounding box: xmin: -84.31423 ymin: 33.88525 xmax: -75.45997 ymax: 36.58003 Geodetic CRS: NAD27 First 10 features: AREA PERIMETER CNTY_ CNTY_ID FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 1 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 2 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 3 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 4 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 5 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 6 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 7 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 8 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 9 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 10 0.114 1.442 1825 1825 37009 5 1091 1 10 1364 SID79 NWBIR79 geometry 1 0 19 POINT (-81.65795 36.58003) 2 0 19 POINT (-81.6387 36.58003) 3 0 19 POINT (-81.61945 36.58003) 4 0 19 POINT (-81.6002 36.58003) 5 0 19 POINT (-81.58096 36.58003) 6 0 19 POINT (-81.56171 36.58003) 7 0 19 POINT (-81.54246 36.58003) 8 0 19 POINT (-81.52321 36.58003) 9 0 19 POINT (-81.50396 36.58003) 10 0 19 POINT (-81.48471 36.58003) > > r = gdal_read(tif) > gt = c(0,1,0,0,0,1) > gdal_inv_geotransform(gt) [1] 0 1 0 0 0 1 > rc = expand.grid(x=1:3, y = 1:3) > #(xy = xy_from_colrow(rc, gt)) > #xy_from_colrow(xy, gt, inverse = TRUE) > crs <- gdal_crs(tif) > > try(gdal_metadata("foo")) [1] NA > gdal_metadata(tif) [1] "AREA_OR_POINT=Point" > > proc.time() user system elapsed 2.53 0.26 2.79