Skip to content

Instantly share code, notes, and snippets.

View mdsumner's full-sized avatar

Michael Sumner mdsumner

  • Integrated Digital East Antarctica, Australian Antarctic Division
  • Hobart, Australia
View GitHub Profile

It will be different depending on whether s2 is enabled or not. To be sure I would use a map projection (local laea: Lambert Azimuthal Equal Area) (for each individual feature, if their coverage is "large enough", fine for something less than continent-size at a layer level) or use s2 (spherical approximating cells) or an ellipsoidal buffering wrapper more directly. You are supposed to blat the crs metadata if you want buffers in longlat for longlat data in sf, but it won't give real geographic distances for many projections that don't support that) and note that terra is very different - it will always use ellipsoidal methods no matter what crs is in use - but has its own vagaries here). It's largely why I don't use these tools for this kind of question, I like functions that do exactly what they say. Doing metrics in geographic space is what map projections exist for and recent attempts to hide this problem from users injects much uncertainty and new problems that are hard to track.

the units on the


## remotes::install_github("USDAForestService/gdalraster") ## since 2024-05-22
library(gdalraster)
dsn <- "NetCDF:/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202209/oisst-avhrr-v02r01.20220925.nc:sst"
translate(dsn, tif <- tempfile(fileext = ".tif", 
                               tmpdir = "/vsimem"), 
          cl_arg = c("-a_srs", "EPSG:4326", "-of", "GTiff", "-co", "COMPRESS=ZSTD"))
aa <- new(VSIFile, tif)
bytes &lt;- aa$ingest(-1)
dsn <- "/vsizip//vsicurl/https://dapds00.nci.org.au/thredds/fileServer/fj7/Copernicus/Sentinel-2/MSI/L1C/2024/2024-05/50S070E-55S075E/S2A_MSIL1C_20240509T044331_N0510_R047_T43FDB_20240509T064234.zip/S2A_MSIL1C_20240509T044331_N0510_R047_T43FDB_20240509T064234.SAFE/GRANULE/L1C_T43FDB_A046376_20240509T044331/IMG_DATA/T43FDB_20240509T044331_TCI.jp2"
system(sprintf("gdalinfo %s", dsn))
Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: /vsizip//vsicurl/https://dapds00.nci.org.au/thredds/fileServer/fj7/Copernicus/Sentinel-2/MSI/L1C/2024/2024-05/50S070E-55S075E/S2A_MSIL1C_20240509T044331_N0510_R047_T43FDB_20240509T064234.zip/S2A_MSIL1C_20240509T044331_N0510_R047_T43FDB_20240509T064234.SAFE/GRANULE/L1C_T43FDB_A046376_20240509T044331/IMG_DATA/T43FDB_20240509T044331_TCI.jp2
Size is 10980, 10980

Python 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0] on linux Type "help", "copyright", "credits" or "license" for more information.

from osgeo import gdal Traceback (most recent call last): File "", line 1, in ModuleNotFoundError: No module named 'osgeo'

set.seed(1)

url = "https://naciscdn.org/naturalearth/10m/raster/SR_LR.zip/SR_LR.tif"
url = paste0("/vsizip/vsicurl/", url)
#url <- "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif"

n = 44000
pts = cbind(x = runif(n, -180, 180), y = runif(n, -90, 90))
## whereever SGAT code does this: 


library(maptools)
data("wrld_simpl")
## Cells per degree
n <- 1
grid <- raster(nrows=n*diff(ylim),ncols=n*diff(xlim),
  xmn=xlim[1],xmx=xlim[2],ymn=ylim[1],ymx=ylim[2],
gdalinfo "/vsicurl/https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5" -sd 26
Driver: HDF5Image/HDF5 Dataset
Files: /vsicurl/https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
Size is 608, 896
Metadata:
  Conventions=CF-1.6
  HDFEOS_INFORMATION_HDFEOSVersion=HDFEOS_5.1.15
  history=This version of the Sea Ice processing code contains updates provided by the science team on September 16, 2019. For details on these updates, see the release notes provided in the DAP.
  institution=NASA's AMSR Science Investigator-led Processing System (SIPS)