Skip to content

Instantly share code, notes, and snippets.

View mdsumner's full-sized avatar

Michael Sumner mdsumner

501 followers · 1.3k following
  • Integrated Digital East Antarctica, Australian Antarctic Division
  • Hobart, Australia
View GitHub Profile
@mdsumner
mdsumner / timing_test_california_wa.md
Last active April 5, 2024 02:30 
dsn <- "/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/L/AC/2023/1/S2A_T01LAC_20230103T221938_L2A/TCI.tif"
library(vapour)
system.time(d <- gdal_raster_data(dsn, bands = 1:3, target_dim = c(1373, 1373)))

0.5 - 2 seconds in California

10-15 seconds in Western Australia or Tasmania

@mdsumner
mdsumner / GADM_geoparquet.md
Last active April 3, 2024 10:26

phwoar, this is a bit exciting, how about convert all those GADM files to one (Geo)Parquet and open it lazily:

url <- "https://geodata.ucdavis.edu/gadm/gadm4.1/gpkg"
f <- readLines(url)
files <- gsub("\"", "", na.omit(stringr::str_extract(f, "gadm.*\\.gpkg\"")))
gpkgs <- file.path(url, files)     


dir.create("raw")
@mdsumner
mdsumner / filelist.R
Created April 3, 2024 05:15
c("/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/L/AH/2023/1/S2B_T01LAH_20230111T222752_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/L/AJ/2023/12/S2B_T01LAJ_20231217T222755_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/K/AB/2023/12/S2A_T01KAB_20231229T221935_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/K/AA/2023/12/S2A_T01KAA_20231229T221935_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/K/AV/2023/12/S2A_T01KAV_20231229T221935_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/K/CB/2023/12/S2B_T01KCB_20231231T220917_L2A/TCI.tif",
"/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/1/K/BB/2023/12/S2B_T01KBB_202
@mdsumner
mdsumner / area_and_crs.md
Last active April 2, 2024 22:59 
    library(terra)
#> terra 1.7.71
set.seed(2)
query <- buffer(project(vect(geosphere::randomCoordinates(1), crs = "EPSG:4326"),   "EPSG:3857"), 
                c(1e6, 4e6), capstyle = "square")

maps::map()
plot(project(query, "EPSG:4326"), add = TRUE, border= "firebrick", lwd = 4)
@mdsumner
mdsumner / kerchunk_serialize_xarray.md
Last active April 1, 2024 20:05 
import fsspec
import kerchunk.hdf

fs = fsspec.filesystem('https')

urls = ["https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240119.nc",
        "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240120.nc",
        "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240121.nc",
        "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202401/oisst-avhrr-v02r01.20240122.nc",
@mdsumner
mdsumner / tiled_window_read_vector_extract.md
Last active March 30, 2024 08:31 
dsn <- "/vsicurl/https://data.source.coop/vizzuality/hfp-100/hfp_2021_100m_v1-2_cog.tif"
dsn <- basename(dsn)
info <- vapour::vapour_raster_info(dsn)

library(grout)
g <- grout(info$dimension, info$extent, blocksize = c(2048, 2048))

index <- tile_index(g)
@mdsumner
mdsumner / spherical.nc.md
Last active March 22, 2024 08:21

This NetCDF has a spherical CRS longlat crs in "crs_wkt", but the NetCDF driver hardcodes WGS84 as the geolocation SRS:

gdalinfo /vsicurl/https://github.com/rgdal-dev/rasterwise/raw/master/data-raw/spherical.nc

...
Geolocation:
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  X
...
@mdsumner
mdsumner / colorize.md
Created March 20, 2024 10:59 
#dsn <- dsn::vrtcon(sds::gebco(), ovr = 4)
dsn <- "vrt:///vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif?ovr=4"

scaled.vrt <- "scaled.vrt"
system(glue::glue("gdal_translate {dsn} scaled.vrt -scale -ot Byte"))

## here I had to edit the pal.vrt to have "scaled.vrt" rather than "mem" in the source (...)
system(glue::glue("gdalattachpct.py pal scaled.vrt pal.vrt -of VRT"))
@mdsumner
mdsumner / ghrsst_R.md
Last active March 19, 2024 23:38

Note that the authorization mechanism requires GDAL 3.6 or later, find your version in terra with terra::gdal().

library(terra)
url <- "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20240318090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"

## make a GDAL dsn for an online file
dsn <- sprintf("/vsicurl/%s", url)
@mdsumner
mdsumner / overviews_.md
Created March 19, 2024 23:17 
gdal_create -outsize 36000 17999 -ot Int8 -co SPARSE_OK=YES -a_srs EPSG:4326 -a_ullr 0 17999 0 36000 weird.tif
gdal_translate weird.tif cog.tif  -of COG
gdalinfo cog.tif | grep Overviews
#  Overviews: 18000x8999, 9000x4499, 4500x2249, 2250x1124, 1125x562, 562x281, 281x14


import xarray