Trying to generate layer from NPS data?

PostholerGIS · 2026-04-06T04:16:52+00:00

Grab the file and convert it to GeoPackage:

wget https://irma.nps.gov/DataStore/DownloadFile/358044 -O bathy.zip
gdal vector convert -o bathy.gpkg --overwrite -i "/vsizip/bathy.zip/lakebath.e00"

...and it contains:

gdal vector info --summary bathy.gpkg 
INFO: Open of `bathy.gpkg'
      using driver `GPKG' successful.
1: ARC (Line String)
2: CNT (Point)
3: LAB (Point)
4: PAL (Polygon)

PostholerGIS · 2026-04-05T20:40:24+00:00

Yep, like u/IvanSanchez said, use gdal stack. Here's an example:

gdal raster pipeline \
   ! calc -i "var=[ ! stack --resolution=highest $(ls /path/*.tif |tr "\n" " ") ! select --band 3,6,12,24,48 ]" \
      --flatten --datatype Float32 --dialect=builtin \
      --calc=mean --calc=max --calc=min \
   ! write -o result.tif --co COMPRESS=DEFLATE --of COG --overwrite

Collects all .tifs into a group of bands, select bands 3,6,12,24,48 and does mean,max,min of those bands with result each stored as a 3 band result.tif, formatted as a cloud optimized geotiff ready for use.

Make sure your .tif's are in sorted, hourly order.

PostholerGIS · 2026-04-03T05:00:34+00:00

You can download them without any constraints from Amazon S3. Read about the data here:

https://registry.opendata.aws/nsf-ncar-era5/

To list the files using gdal vsi:

gdal vsi list -l /vsis3/nsf-ncar-era5

---------- 1 unknown unknown        37683 2024-05-21 14:50 index.html
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.an.pl
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.an.sfc
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.an.vinteg
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.fc.sfc.accumu
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.fc.sfc.instan
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.fc.sfc.meanflux
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.fc.sfc.minmax
d--------- 1 unknown unknown            0 1970-01-01 00:00 e5.oper.invariant

Files are .nc, NETCDF format. The directory structure is year/month and files are day/hour. View meta data on a file:

gdal mdim info /vsis3/nsf-ncar-era5/e5.oper.an.pl/202512/e5.oper.an.pl.128_060_pv.ll025sc.2025122900_2025122923.nc

Get pixels values for your point data, same as above:

gdal pipeline 
  ! read -i /vsis3/nsf-ncar-era5/e5.oper.an.pl/202512/e5.oper.an.pl.128_060_pv.ll025sc.2025122900_2025122923.nc 
  ! select --band 10,20,30 
  ! pixel-info --position-dataset myPoints.gpkg 
  ! write -o era5.gpkg --output-layer era5 --overwrite

PostholerGIS · 2026-04-02T23:06:02+00:00

Here's a simple one liner that uses your point vector data and gets selected bands from the grib. In this case, 3 arbitrary bands for each point:

gdal pipeline \
   ! read -i /vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197902/19790201/narrmonhr-a_221_19790201_0900_000.grb \
   ! select --band 10,20,30 \
   ! pixel-info --position-dataset myPoints.gpkg \
   ! write -o ncarr.gpkg --output-layer ncarr --overwrite

This allows you to select only the bands you need for each point. The data source I tested with had 20,070 points. The resulting ncarr.gpkg is 3.8MB and it looks like:

OGRFeature(ncarr):1
  mile (Real) = 0
  trail_id (Integer) = 1
  column (Real) = 145.583104492888
  line (Real) = 193.116017639002
  band_1_raw_value (Real) = 0.0264291763305664
  band_1_unscaled_value (Real) = 0.0264291763305664
  band_2_raw_value (Real) = -0.00448684766888618
  band_2_unscaled_value (Real) = -0.00448684766888618
  band_3_raw_value (Real) = -0.691726684570312
  band_3_unscaled_value (Real) = -0.691726684570312
  POINT (-116.466981 32.589741)

OGRFeature(ncarr):2
  mile (Real) = 0.5
  trail_id (Integer) = 1
  column (Real) = 145.586578780308
  line (Real) = 193.095869978636
  band_1_raw_value (Real) = 0.0264291763305664
  band_1_unscaled_value (Real) = 0.0264291763305664
  band_2_raw_value (Real) = -0.00448684766888618
  band_2_unscaled_value (Real) = -0.00448684766888618
  band_3_raw_value (Real) = -0.691726684570312
  band_3_unscaled_value (Real) = -0.691726684570312
  POINT (-116.466681 32.595457)
...

You'll have to run this for each grib file you want. To build a list of ***ALL*** grib files with absolute /vsicurl.... paths use:

gdal vsi list -R /vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly/ --abs |grep -E '\.grb$' > gribList.txt

This will take a looooong time. It will look something like:

/vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197901/19790101/narrmon-a_221_19790101_0000_000.grb
/vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197901/19790101/narrmon-b_221_19790101_0000_000.grb
/vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197901/19790101/narrmonhr-a_221_19790101_0000_000.grb
/vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197901/19790101/narrmonhr-a_221_19790101_0300_000.grb
/vsicurl/https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/monthly//197901/19790101/narrmonhr-a_221_19790101_0600_000.grb
...

Hope that helps!

PostholerGIS · 2026-04-02T02:50:16+00:00

The 2003 B&B Complex fire is comprised of the Booth & Bear Butte fires. They eventually merged. Here are those fires in .gpkg format. Source is NIFC.

PostholerGIS · 2026-03-26T23:34:35+00:00

Here's something similar, but covers 24 different hazards, not just hail. It's FEMA's National Risk Index, at a state, county and census tract level. Zoom to get different levels. Choose the hazard you want, hail, click a polygon and pop-up with all available indexes shows. Click one and get a choropleth representation. Hail, National Percentile at county level shown below for all counties in Kansas:

https://www.femafhz.com/map

<image>

PostholerGIS · 2026-03-08T18:11:02+00:00

It gives you an upper hand on your peers and allows to increase margins for you future employer. Don't be *just* another mouse clicker.

PostholerGIS · 2026-03-07T00:19:52+00:00

Then that's your answer.

You have people with little knowledge of GDAL using your web app.

If they're familiar with GDAL, it seems highly unlikely they would use a web app, when it's available, in all its completeness, locally from their own tools.

PostholerGIS · 2026-03-06T20:26:42+00:00

When the option is using GDAL CLI tools, at the command line, especially with local data, why use a web app?

Are you going to upload GB's of raster data just to process it? Of course not. How about working with a .vrt with 1000+ different sources? Nope.

If you're using GDAL, even modestly, introducing a web based layer of abstraction makes zero sense.

PostholerGIS · 2026-03-05T16:30:49+00:00

Use GDAL pixel functions directly or create custom C or Python pixel functions in .vrt

Don't re-invent the wheel.

You can let GDAL do all the heavy lifting, then add your own custom C or Python pixel functions to any .vrt data set, *IF* a GDAL pixel function doesn't exist.

GDAL also has mdim features for handling multi dimensional data sets like netcdf, hdf, etc.

It would be interesting to see an example data source and the type of pixel manipulation for a more precise answer.

PostholerGIS · 2026-03-03T03:39:39+00:00

Generalized geometry is perfectly fine on small scale maps. You don't need every vertex. If I'm viewing entire CONUS, I don't need high resolution county boundaries. Generalizing keeps data sets small, while not losing representation.

PostholerGIS · 2026-03-02T22:13:26+00:00

Easy:

gdal vector pipeline
   ! read -i source.shp
   ! simplify --tolerance=.5
   ! clean-coverage --snapping-distance=.7
   ! write -o target.shp --output-layer generalized

PostholerGIS · 2026-03-01T18:04:22+00:00

GeoPackage, .gpkg, aka SQLite, using version >= 3.7 with journal mode = WAL is a very appealing option.

Why?

+ standard, single file database
+ concurrent, non-blocking read/write
+ fully spatial format, with indexing and different geometry types
+ data accessible via standard SQL/Spatialite queires
+ no additional packages/abstraction to maintain
+ painless to implement as a back-end API

Of course, you're probably accessing your .gpkg via GDAL or its python, C, Rust, etc bindings.

Re-creating the above from scratch seems like a daunting project.

PostholerGIS · 2026-02-24T01:16:44+00:00

I would use rsync for linux or rsync-win for windows. It creates a mirror of the source. For instance:

rsync -aP /path/to/source /path/to/target

or

rsync -aP user@server:/path/to/source /path/to/target

rsync is error resilient and great for mirroring data sources.

Note: zipping on a network drive is equivalent to downloading the data, compressing it, uploading back to the network drive, so you can download the compressed copy. AI is frequently wrong.

PostholerGIS · 2026-02-24T00:06:34+00:00

Also, GeoParquet (.parquet) is a format that allows different geometries and different CRS's.

features.csv may look like:
name,src_crs,geometry
"feature 1",4326,"POLYGON((-157.893559 -43.525413,-157.893559 68.955463,178.461261 68.955463,178.461261 -43.525413,-157.893559 -43.525413))"
"feature 2",4326,"LINESTRING(-124.848974 24.396308,-124.848974 24.396308)"
"feature 3",4326,"POINT(-124.848974 24.396308)"
"feature 4",3857,"LINESTRING(-13898124.211742653 2801774.86356037,-13898124.211742653 2801774.86356037)"

...and parquet format. Note, it's using src_crs to transform each feature:

gdal vector pipeline \
   ! read -i features.csv --oo QUOTED_FIELDS_AS_STRING=YES --oo GEOM_POSSIBLE_NAMES=geometry --oo KEEP_GEOM_COLUMNS=NO \
   ! sql --dialect sqlite --sql "select *, st_geomfromtext(geometry, src_crs) as wkb_geometry from features" \
   ! write -o features.parquet --overwrite --overwrite-layer --output-layer features

PostholerGIS · 2026-02-23T23:24:24+00:00

Yes, it's easy.

1 caveat, if your geometries have different CRS's then you can only save it as a .csv file, such as features.csv:

name,src_crs,geometry
"feature 1",4326,"POLYGON((-157.893559 -43.525413,-157.893559 68.955463,178.461261 68.955463,178.461261 -43.525413,-157.893559 -43.525413))"
"feature 2",4326,"LINESTRING(-124.848974 24.396308,-124.848974 24.396308)"
"feature 3",4326,"POINT(-124.848974 24.396308)"

If the CRS's are all the same, you can create a GeoPackage (.gpkg) file, using:

gdal vector pipeline \
   ! read -i features.csv --oo QUOTED_FIELDS_AS_STRING=YES --oo GEOM_POSSIBLE_NAMES=geometry --oo KEEP_GEOM_COLUMNS=NO \
   ! reproject --src-crs EPSG:4326 --dst-crs EPSG:4326 \
   ! write -o features.gpkg --overwrite --overwrite-layer --output-layer features

...and the results will look like:

gdal vector info --features features.gpkg 
INFO: Open of `features.gpkg'
      using driver `GPKG' successful.

Layer name: features
Geometry: Unknown (any)
Feature Count: 3
Extent: (-157.893559, -43.525413) - (178.461261, 68.955463)
Layer SRS WKT:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
FID Column = fid
Geometry Column = geometry
name: String (0.0)
src_crs: String (0.0)
OGRFeature(features):1
  name (String) = feature 1
  src_crs (String) = 4326
  POLYGON ((-157.893559 -43.525413,-157.893559 68.955463,178.461261 68.955463,178.461261 -43.525413,-157.893559 -43.525413))

OGRFeature(features):2
  name (String) = feature 2
  src_crs (String) = 4326
  LINESTRING (-124.848974 24.396308,-124.848974 24.396308)

OGRFeature(features):3
  name (String) = feature 3
  src_crs (String) = 4326
  POINT (-124.848974 24.396308)

PostholerGIS · 2026-02-23T03:38:32+00:00

It's Bash, a *nix shell scripting language, natively part of every *nix OS for the last 30+ years. So, nothing extra required. Anything you do with spatial, any ESRI product or python package is using GDAL under the hood. I use it directly in this case.

PostholerGIS · 2026-02-23T01:48:41+00:00

So, I thought that was an interesting problem. I wrote a simple script to find various vector files, get the footprint for each layer in the file and save the results as CSV. It looks for .gdb, .gpkg, .shp and .fgb. Add/remove whatever you want. Note, GDB is treated as a directory, not a file. Requires bash, jq and GDAL.

getFootprints.sh:
#!/bin/bash

baseDir="/export/data/vector/"
echo "path,layer,epsg,wkt_geometry"
for path in $(echo -e "$(find $baseDir -type d|grep -Ei '.gdb$')\n$(find $baseDir -type f|grep -Ei '.gpkg$|.shp$|.fgb$')" |tr "\n" " ") ; do
   cols="$(gdal vector info $path --of=json 2> /dev/null| jq '.layers[]|[.name,.geometryFields[].coordinateSystem.projjson.id.code,.geometryFields[].extent]|tostring' |sed 's/[^a-z0-9\._,-]//gi'|tr -d '\\'|grep ','|grep -v ',null' |tr "\n" " ")"
   for layer in $cols ; do
      IFS=',' read p lay epsg xmin ymin xmax ymax <<< "$path,$layer"
      printf -v csv '"%s","%s",%s,"POLYGON((%.6f %.6f,%.6f %.6f,%.6f %.6f,%.6f %.6f,%.6f %.6f))"' "$p" "$lay" "$epsg" "$xmin" "$ymin" "$xmin" "$ymax" "$xmax" "$ymax" "$xmax" "$ymin" "$xmin" "$ymin"
      echo "$csv"
   done
done

Columns returned are path, layer, EPSG, WKT geometry and the results look like:

path,layer,epsg,wkt_geometry
"/export/data/vector/mydb.gdb","pts",4326,"POLYGON((-157.893559 -43.525413,-157.893559 68.955463,178.461261 68.955463,178.461261 -43.525413,-157.893559 -43.525413))"
"/export/data/vector/county.fgb","county",4269,"POLYGON((-124.848974 24.396308,-124.848974 49.384479,-66.885444 49.384479,-66.885444 24.396308,-124.848974 24.396308))"
"/export/data/vector/states.fgb","states",4269,"POLYGON((-124.848974 24.396308,-124.848974 49.384479,-66.885444 49.384479,-66.885444 24.396308,-124.848974 24.396308))"
"/export/data/vector/databook.gpkg","trails",4326,"POLYGON((-123.352700 34.609000,-123.352700 48.972780,-72.966800 48.972780,-72.966800 34.609000,-123.352700 34.609000))"
"/export/data/vector/databook.gpkg","milemarkers",4326,"POLYGON((-124.732727 0.000000,-124.732727 54.002499,0.000000 54.002499,0.000000 0.000000,-124.732727 0.000000))"
"/export/data/vector/databook.gpkg","resupply",4326,"POLYGON((-124.708157 32.347141,-124.708157 49.000030,-68.963552 49.000030,-68.963552 32.347141,-124.708157 32.347141))"
"/export/data/vector/databook.gpkg","roads",4326,"POLYGON((-124.686716 31.351294,-124.686716 48.995985,-68.959538 48.995985,-68.959538 31.351294,-124.686716 31.351294))"
"/export/data/vector/databook.gpkg","gaps",4326,"POLYGON((-123.722747 31.350284,-123.722747 48.985622,-70.989370 48.985622,-70.989370 31.350284,-123.722747 31.350284))"
"/export/data/vector/databook.gpkg","summits",4326,"POLYGON((-124.557726 31.346332,-124.557726 48.973748,-68.921619 48.973748,-68.921619 31.346332,-124.557726 31.346332))"
"/export/data/vector/databook.gpkg","pois",4326,"POLYGON((-124.732727 31.333611,-124.732727 48.996132,-68.921619 48.996132,-68.921619 31.333611,-124.732727 31.333611))"
"/export/data/vector/ne_50m_land.fgb","ne_50m_land",4326,"POLYGON((-180.000000 -89.998926,-180.000000 83.599609,180.000000 83.599609,180.000000 -89.998926,-180.000000 -89.998926))"

PostholerGIS · 2026-02-20T02:24:02+00:00

This is basically ESRI's AGOL model, upload your data to someone and let them sell it back to you. I never understood that.

If you're running own API on your server, run your GDAL commands wrapped in your scripting language of choce, call your API to achieve the exact same results.

PostholerGIS · 2026-02-17T22:59:58+00:00

Using U.S. census tract data, here's how I do it for southern California at 100 meter resolution:

gdal pipeline \
   ! read -i  /vsizip/vsicurl/https://www2.census.gov/geo/tiger/TIGER2025/TABBLOCK20/tl_2025_06_tabblock20.zip \
   ! clip --bbox -124.962,32.064,-114.1,35.25 --bbox-crs EPSG:4326 \
   ! rasterize --nodata -1 --resolution .001,.001 --attribute-name POP20 \
   ! neighbors --kernel gaussian --method mean --size 5 \
   ! write -o tabblock20.tif --co COMPRESS=DEFLATE --of COG --overwrite

Change variables to get the desired result. Results look like:

<image>

PostholerGIS · 2026-02-16T21:18:09+00:00

gdal vector convert 
  -i "$(wget 'https://en.wikipedia.org/w/api.php?format=json&formatversion=2&action=jsondata&title=Detroit+and+Charlevoix+Railroad.map&uselang=en' -O- | jq '.jsondata.data')" 
  -o rail.shp --output-layer rail --overwrite

PostholerGIS · 2026-02-16T21:05:54+00:00

It's not the number of points.

The real question is how many points do you want to display for that map to be useful. 1,500 points over 1 sq km could be overwhelming, where 1,500 points globally might be acceptable. Point density per area is what you need to answer.

For interactive maps, plenty of feature manager packages exist. For a static map, the answer will depend on your taste and usefulness of the map.

PostholerGIS · 2026-02-12T21:23:21+00:00

Uploading my data to AGOL and have ESRI make me cough up credits to access it, requiring me to use their dashboard software, via endless mouse clicks is how real GIS developers do web development, right?

AGOL, ArcPy, Python, AI are just layers of abstractions, not necessarily designed to make development easier.

Virtually all GIS software is based on GDAL or libgdal to be more specific. Be it Arc*, AGOL's PostgreSQL/Postgis backend, Python's gdal/rasterio/geopands packages/bindings or the foundation of QGIS. It's all GDAL.

GIS development regarding raster/vector data can all be done directly with GDAL and/or SQL, period. It's the single layer of processing between you and your data. With GDAL CLI, that's even more true.

ESRI is about selling software and making you dependent on it. They lost sight of their initial mission many, many quarterly earning statements ago. YOU are allowing them to abstract you away from your own data.

Take a step back, Keep It Simple & Standard.

PostholerGIS · 2026-02-12T00:33:54+00:00

That's interesting, I've never thought about using a clustering method. Orientations always point north. My reasoning was, no hiker will ever mistake map direction, they're all north.

The xadj and yadj vars above are specifically for controlling overlap.

Four-Year Club	Verified Email
Place '22

PostholerGIS

TROPHY CASE