mdsumner

ATLAS/ICESat-2 L3B Daily and Monthly Gridded Sea Ice Freeboard

https://nsidc.org/data/atl20/versions/4

https://n5eil01u.ecs.nsidc.org/ATLAS/ATL20.004/2018.10.14/ATL20-01_20181001010615_00370101_004_01.h5

there’s a problem with the grid_mapping, the ‘crs’ declaration is non-compliant and prevents downstream tools picking up the projection and georeferenced metadata.

ncdump -h ATL20-02_20181001010615_00370101_004_01.h5 | grep grid_mapping

mdsumner

## 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],

mdsumner

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)

mdsumner

user guide: https://nsidc.org/sites/default/files/au_si12-v001-userguide_1.pdf (lists the extent for the grid, it's the same as the 25km ones)

https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/

Can read this with your earthdata token set up

header <- "Authorization: Bearer <token>"
u <- "/vsicurl/https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5"
Sys.setenv("GDAL_HTTP_HEADERS" = header)

mdsumner

See

> ## Number of chunks and runs per chunk
> chunks <- 250
> runs <- 40
> plan(multisession, workers= availableCores()-1)
Error: Cannot create 127 parallel PSOCK nodes. Each node needs one connection, but there are only 124 connections left out of the maximum 128 available on this R installation

mdsumner

GeoBox from a DSN

def extent(dimension, transform): 
  return [transform[0], transform[0] + dimension[0] * transform[1], 
  transform[3] + dimension[1] * transform[5], transform[3]]
  
  
def bbox(dimension, transform): 

mdsumner

library(duckdbfs)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 

mdsumner

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.9.0dev-cb4d30f56d, PROJ 9.3.1; sf_use_s2() is
#> TRUE
library(raster)
#> Loading required package: sp


## threshold is roughly side length in metres of the roi, i.e. sqrt(area)
dem_source <- function(threshold = 1) {

mdsumner

I don't entirely trust the timings here, but it's faster (maybe 2x faster to plot).

  dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"

library(gdalraster)
#> GDAL 3.9.0dev-cb4d30f56d, released 2024/04/15 (debug build), GEOS 3.12.1, PROJ 9.3.1


ds_int <- new(GDALRaster, sprintf("vrt://%s?ot=Int32", dsn))

mdsumner

from osgeo import gdal
gdal.UseExceptions()



## we do something with vrt
dsn = "vrt:///vsicurl/https://github.com/cran/rgl/raw/master/inst/textures/worldsmall.png?a_ullr=-180,90,180,-90&expand=rgb"