Flynn et al. (2026) created Stereo DEMS of the 2022 Mauna Loa eruption and differenced them against pre-eruption USGS LiDAR DEMs to quatify lava flow thickness. This notebook contains code to convert the original USGS LiDAR DEMs from their native CRS (EPSG: 6318) to UTM WGS2139, which is the CRS of the post-eruption stereo DEMs. Please refer to the publication that describes the methods and results in detail!
Method Overview:
Geoid to Ellipsoid: Reproject original EPSG:6635 [NAD83(PA11) / UTM zone 5N] to EPSG:6321 [NAD83(PA11)]
Datum Shift: Intermediate EPSG:6321 [NAD83(PA11)] to UTM_Zone5_G2139 [WGS84 G2139 with custom WKT]
Mosaic: Multiple input DEMs to match post-eruption stereo DEM grid for differencing
Background¶
Input DEMs were obtained from USGS 3DEP Lidar via the OpenTopgraphy API:
https://
According to the collection metadata, this LiDAR data is referenced to GEOID12B. Shifting to ellipsoid heights is done empirically with a vertical shift grid available at https://
import numpy as np
import rioxarray
import pyprojStage data¶
This notebook reprojects several large DEMs (~700Mb) stored as GeoTiffs. In the following code we download a folder of all DEMs (~2 GB) to process with GDAL
%%bash
# This will take ~2 min on typical WiFi connections
cd /tmp
wget -nc --no-verbose https://github.com/uw-cryo/3D_CRS_Transformation_Resources/releases/download/v0.1/flynn_jvgr2026_3dep.zip
unzip -n flynn_jvgr2026_3dep.zip
ls -lh flynn_jvgr2026_3dep/Archive: flynn_jvgr2026_3dep.zip
total 8967048
-rw-r--r--@ 1 scotthenderson wheel 722M Mar 23 09:34 Kokoolau_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 722M Mar 23 09:32 MaunaLoa_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 722M Mar 23 09:35 PuuKoli_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 715M Mar 23 09:36 PuuOo_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 722M Mar 23 09:36 PuuUlaula_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 722M Mar 23 09:32 SulphurCone_3DEP.tif
-rw-r--r--@ 1 scotthenderson wheel 1.5K Mar 23 14:01 UTM_Zone5_G2139.wkt
Review original CRS¶
Unfortunately, the original metadata did not record the vertical datum, but reading the 3DEP documentation data over Hawaii references the GEOID12B model. Read more at https://
# Ensure we're working in the data folder
%cd /tmp/flynn_jvgr2026_3dep/private/tmp/flynn_jvgr2026_3dep
ml_lidar_fn = 'MaunaLoa_3DEP.tif'
da = rioxarray.open_rasterio(ml_lidar_fn,masked=True)
crs = da.rio.crs# Original 2D CRS
pyproj.CRS(crs)<Projected CRS: EPSG:6635>
Name: NAD83(PA11) / UTM zone 5N
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: UTM zone 5N
- method: Transverse Mercator
Datum: NAD83 (National Spatial Reference System PA11)
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich# Our target: Lon,Lat Ellipsoid height
pyproj.CRS('EPSG:6321')<Geographic 3D CRS: EPSG:6321>
Name: NAD83(PA11)
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: American Samoa, Marshall Islands, United States (USA) - Hawaii, United States minor outlying islands; onshore and offshore.
- bounds: (157.47, -17.56, -151.27, 31.8)
Datum: NAD83 (National Spatial Reference System PA11)
- Ellipsoid: GRS 1980
- Prime Meridian: GreenwichWarp with GDAL¶
Step 1: Convert to ellipsoid height¶
force a specific geoid shift grid (us_noaa_g2012bh0.tif) by specifying PROJ pipeline:
+proj=pipeline
+step +inv +proj=utm +zone=5 +ellps=GRS80
+step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1
+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
+step +proj=axisswap +order=2,1%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=MaunaLoa_3DEP.tif
OUTPUT=MaunaLoa_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13905P x 13862L.
Using internal nodata values (e.g. 3.4e+38) for image MaunaLoa_3DEP.tif.
Copying nodata values from source MaunaLoa_3DEP.tif to destination MaunaLoa_3DEP_6321.tif.
Processing MaunaLoa_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:07.
Step 2: Reproject to WGS84 G2139 UTM¶
EPSG codes for UTM do not reference specific WGS84 realizations, so we use a custom CRS definition to ensure proper reprojection and recording of the vertical datum.
Once we have data in EPSG:6321 (lat/lon, heights above ellipsoid), we can go into our UTM definition with properly specified datum realization (G2139) with very high accuracy (0.01 m).
# This uses a Helmert transform to go from NAD83 -> WGS84
!projinfo -o proj -s EPSG:6321 -t "$(cat UTM_Zone5_G2139.wkt)" --bbox -155.5727033,19.4101513,-155.4696472,19.7729149 --spatial-test intersects --hide-ballpark
Candidate operations found: 1
-------------------------------------
Operation No. 1:
unknown id, Conversion from NAD83(PA11) (geog3D) to NAD83(PA11) (geocentric) + Inverse of ITRF2014 to NAD83(PA11) (1) + Inverse of WGS 84 (G2139) to ITRF2014 (1) + Conversion from WGS 84 (G2139) (geocentric) to WGS 84 (G2139) (geog3D) + UTM zone 5N, 0.01 m, American Samoa, Marshall Islands, United States (USA) - Hawaii, United States minor outlying islands; onshore and offshore., time-dependent operation
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
+step +proj=cart +ellps=GRS80
+step +inv +proj=helmert +x=0.9109 +y=-2.0129 +z=-0.5863 +rx=0.022749
+ry=0.02656 +rz=-0.025706 +s=0.00212 +dx=0.0001 +dy=0.0001 +dz=-0.0019
+drx=-0.000384 +dry=0.001007 +drz=-0.002186 +ds=0.00011 +t_epoch=2010
+convention=coordinate_frame
+step +inv +proj=cart +ellps=WGS84
+step +proj=utm +zone=5 +ellps=WGS84
!gdalwarp -wo NUM_THREADS=ALL_CPUS -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt MaunaLoa_3DEP_6321.tif MaunaLoa_3DEP_32605_G2139.tifCreating output file that is 13756P x 14435L.
Using internal nodata values (e.g. 3.4e+38) for image MaunaLoa_3DEP_6321.tif.
Copying nodata values from source MaunaLoa_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/MaunaLoa_3DEP_32605_G2139.tif_21448_1.
0...10]9;4;1;150...10.. - estimated remaining time: 00:00:2]9;4;1;170...10... - estimated remaining time: 00:00:2]9;4;1;200...10...20 - estimated remaining time: 00:00:1]9;4;1;220...10...20. - estimated remaining time: 00:00:2]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:1]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:1]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:1]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:27.
Sulphur Cone DEM Conversion¶
%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=SulphurCone_3DEP.tif
OUTPUT=SulphurCone_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13927P x 13883L.
Using internal nodata values (e.g. 3.4e+38) for image SulphurCone_3DEP.tif.
Copying nodata values from source SulphurCone_3DEP.tif to destination SulphurCone_3DEP_6321.tif.
Processing SulphurCone_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:06.
!gdalwarp -wo NUM_THREADS=ALL_CPUS -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt SulphurCone_3DEP_6321.tif SulphurCone_3DEP_32605_G2139.tifCreating output file that is 13789P x 14468L.
Using internal nodata values (e.g. 3.4e+38) for image SulphurCone_3DEP_6321.tif.
Copying nodata values from source SulphurCone_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/SulphurCone_3DEP_32605_G2139.tif_24545_1.
0...10]9;4;1;150...10.. - estimated remaining time: 00:00:2]9;4;1;170...10... - estimated remaining time: 00:00:2]9;4;1;200...10...20 - estimated remaining time: 00:00:2]9;4;1;220...10...20. - estimated remaining time: 00:00:1]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:1]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:1]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:1]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:27.
Kokoolau DEM Conversion¶
%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=Kokoolau_3DEP.tif
OUTPUT=Kokoolau_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13903P x 13860L.
Using internal nodata values (e.g. 3.4e+38) for image Kokoolau_3DEP.tif.
Copying nodata values from source Kokoolau_3DEP.tif to destination Kokoolau_3DEP_6321.tif.
Processing Kokoolau_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:06.
!gdalwarp -wo NUM_THREADS=ALL_CPUS -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt Kokoolau_3DEP_6321.tif Kokoolau_3DEP_32605_G2139.tifCreating output file that is 13750P x 14439L.
Using internal nodata values (e.g. 3.4e+38) for image Kokoolau_3DEP_6321.tif.
Copying nodata values from source Kokoolau_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/Kokoolau_3DEP_32605_G2139.tif_27606_1.
0...10]9;4;1;150...10.. - estimated remaining time: 00:00:2]9;4;1;170...10... - estimated remaining time: 00:00:2]9;4;1;200...10...20 - estimated remaining time: 00:00:2]9;4;1;220...10...20. - estimated remaining time: 00:00:2]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:1]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:1]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:0]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:25.
PuuKoli DEM Conversion¶
%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=PuuKoli_3DEP.tif
OUTPUT=PuuKoli_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13900P x 13857L.
Using internal nodata values (e.g. 3.4e+38) for image PuuKoli_3DEP.tif.
Copying nodata values from source PuuKoli_3DEP.tif to destination PuuKoli_3DEP_6321.tif.
Processing PuuKoli_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:06.
!gdalwarp -wo NUM_THREADS=ALL_CPUS -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt PuuKoli_3DEP_6321.tif PuuKoli_3DEP_32605_G2139.tifCreating output file that is 13744P x 14443L.
Using internal nodata values (e.g. 3.4e+38) for image PuuKoli_3DEP_6321.tif.
Copying nodata values from source PuuKoli_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/PuuKoli_3DEP_32605_G2139.tif_30477_1.
0...10]9;4;1;170...10... - estimated remaining time: 00:00:2]9;4;1;200...10...20 - estimated remaining time: 00:00:1]9;4;1;220...10...20. - estimated remaining time: 00:00:1]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:0]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:0]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:0]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:25.
PuuOo DEM Conversion¶
%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=PuuOo_3DEP.tif
OUTPUT=PuuOo_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13879P x 13837L.
Using internal nodata values (e.g. 3.4e+38) for image PuuOo_3DEP.tif.
Copying nodata values from source PuuOo_3DEP.tif to destination PuuOo_3DEP_6321.tif.
Processing PuuOo_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:06.
!gdalwarp -wo NUM_THREADS=ALL_CPUS -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt PuuOo_3DEP_6321.tif PuuOo_3DEP_32605_G2139.tifCreating output file that is 13712P x 14412L.
Using internal nodata values (e.g. 3.4e+38) for image PuuOo_3DEP_6321.tif.
Copying nodata values from source PuuOo_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/PuuOo_3DEP_32605_G2139.tif_33369_1.
0...10]9;4;1;150...10.. - estimated remaining time: 00:00:2]9;4;1;170...10... - estimated remaining time: 00:00:2]9;4;1;200...10...20 - estimated remaining time: 00:00:1]9;4;1;220...10...20. - estimated remaining time: 00:00:2]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:1]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:0]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:0]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:25.
PuuUlaula DEM Conversion¶
%%bash
PROJ_PIPELINE='+proj=pipeline +step +inv +proj=utm +zone=5 +ellps=GRS80 +step +proj=vgridshift +grids=us_noaa_g2012bh0.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1'
INPUT=PuuUlaula_3DEP.tif
OUTPUT=PuuUlaula_3DEP_6321.tif
CPL_DEBUG=OFF PROJ_DEBUG=2 \
PROJ_NETWORK=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -overwrite -wm 500 -r bilinear \
-wo NUM_THREADS=ALL_CPUS \
-co TILED=YES \
-co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
-s_srs EPSG:6635+5703 -t_srs EPSG:6321 \
-ct "${PROJ_PIPELINE}" \
${INPUT} ${OUTPUT}Creating output file that is 13881P x 13840L.
Using internal nodata values (e.g. 3.4e+38) for image PuuUlaula_3DEP.tif.
Copying nodata values from source PuuUlaula_3DEP.tif to destination PuuUlaula_3DEP_6321.tif.
Processing PuuUlaula_3DEP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done in 00:00:05.
!gdalwarp -tap -of COG -r bilinear -tr 1.0 1.0 -s_srs EPSG:6321 -t_srs UTM_Zone5_G2139.wkt PuuUlaula_3DEP_6321.tif PuuUlaula_3DEP_32605_G2139.tifCreating output file that is 13718P x 14409L.
Using internal nodata values (e.g. 3.4e+38) for image PuuUlaula_3DEP_6321.tif.
Copying nodata values from source PuuUlaula_3DEP_6321.tif to destination /var/folders/1v/k85p3x5d10zb9mhqx88lb3240000gn/T/PuuUlaula_3DEP_32605_G2139.tif_36275_1.
0.]9;4;1;70... - estimated remaining time: 00:01:0]9;4;1;100...10 - estimated remaining time: 00:01:0]9;4;1;120...10. - estimated remaining time: 00:01:0]9;4;1;150...10.. - estimated remaining time: 00:01:0]9;4;1;170...10... - estimated remaining time: 00:01:0]9;4;1;200...10...20 - estimated remaining time: 00:01:0]9;4;1;220...10...20. - estimated remaining time: 00:00:5]9;4;1;250...10...20.. - estimated remaining time: 00:00:5]9;4;1;270...10...20... - estimated remaining time: 00:00:5]9;4;1;300...10...20...30 - estimated remaining time: 00:00:5]9;4;1;320...10...20...30. - estimated remaining time: 00:00:5]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:4]9;4;1;370...10...20...30... - estimated remaining time: 00:00:4]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:4]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:4]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:4]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:4]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:3]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:3]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:3]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:3]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:3]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:2]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:2]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:2]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:2]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:2]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:1]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:1]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:1]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:1]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:1]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:01:01.
Mosaic the DEMs¶
We combine all DEMs into ML_pre_event_3DEP.tif using Ames Stereo Pipelines dem_mosaic tool, which you can install with:
pixi global install stereo-pipeline -c conda-forge -c nasa-ames-stereo-pipeline "3.5.0_no_isis"
Full mosaic¶
First we create a mosaic using the full merged extent of all input DEMs
!dem_mosaic -o ML_pre_event_3DEP.tif *32605_G2139.tif --> Setting number of processing threads to: 4
Writing log: ML_pre_event_3DEP.tif-log-dem_mosaic-03-24-1153-42067.txt
Using output no-data value: 3.4e+38
Determining the bounding boxes of the inputs.
--> [*************************************************************] 100%
Mosaic size: 40408 x 42123 pixels.
Number of tiles: 1 x 1 = 1.
Reading the input DEMs.
Loading: Kokoolau_3DEP_32605_G2139.tif
Loading: MaunaLoa_3DEP_32605_G2139.tif
Loading: PuuKoli_3DEP_32605_G2139.tif
Loading: PuuOo_3DEP_32605_G2139.tif
Loading: PuuUlaula_3DEP_32605_G2139.tif
Loading: SulphurCone_3DEP_32605_G2139.tif
Writing: ML_pre_event_3DEP.tif
--> [*************************************************************] 100%
Re-writing with blocks of size: 256 x 256.
--> [*************************************************************] 99%
Number of valid (not no-data) pixels written: 1051634157.
Cropped mosaic to match post-even stereo DEM¶
The following functions are used to crop and match the grid of the post lava flow stereo DEMs so that we can take a difference to extract lava flow volume. For which we obtain bounds: ‘217647.0 2145461.0 244140.0 2193138.0’
def nearest_floor(x, a):
"""
Round down to the nearest smaller multiple of a.
From https://github.com/uw-cryo/EarthLab_AirQuality_UAV/blob/main/notebooks/EarthLab_AQ_lidar_download_processing_function.ipynb
"""
return np.floor(x / a) * a
def nearest_ceil(x, a):
"""
Round down to the nearest larger multiple of a.
From https://github.com/uw-cryo/EarthLab_AirQuality_UAV/blob/main/notebooks/EarthLab_AQ_lidar_download_processing_function.ipynb
"""
return np.ceil(x / a) * a
def tap_bounds(site_bounds, res):
"""
Return target aligned pixel bounds for a given site bounds and resolution.
From https://github.com/uw-cryo/EarthLab_AirQuality_UAV/blob/main/notebooks/EarthLab_AQ_lidar_download_processing_function.ipynb
"""
return ([nearest_floor(site_bounds[0], res), nearest_floor(site_bounds[1], res), \
nearest_ceil(site_bounds[2], res), nearest_ceil(site_bounds[3], res)])
# stereo_dsm = rioxarray.open_rasterio('20230702_20230801_wt_avg_mos.tif')
# stereo_dsm.rio.crs
# bounds = stereo_dsm.rio.bounds()
# bounds_buffered = [bounds[0]-10,bounds[1]-10,bounds[2]+10,bounds[3]+10] #buffer bounds by 10 m
# tapped_bounds = tap_bounds(bounds_buffered,res=1.0)
# str_bounds = [str(x) for x in tapped_bounds]
# t_projwin = ' '.join(str_bounds)!dem_mosaic --t_projwin 217647.0 2145461.0 244140.0 2193138.0 -o ML_pre_event_3DEP_stereo_bounds.tif *32605_G2139.tif --> Setting number of processing threads to: 4
Writing log: ML_pre_event_3DEP_stereo_bounds.tif-log-dem_mosaic-03-24-1207-20514.txt
Using output no-data value: 3.4e+38
Determining the bounding boxes of the inputs.
--> [*************************************************************] 100%
Mosaic size: 26494 x 40349 pixels.
Number of tiles: 1 x 1 = 1.
Reading the input DEMs.
Loading: Kokoolau_3DEP_32605_G2139.tif
Loading: MaunaLoa_3DEP_32605_G2139.tif
Loading: PuuKoli_3DEP_32605_G2139.tif
Loading: PuuOo_3DEP_32605_G2139.tif
Loading: PuuUlaula_3DEP_32605_G2139.tif
Loading: SulphurCone_3DEP_32605_G2139.tif
Writing: ML_pre_event_3DEP_stereo_bounds.tif
--> [*************************************************************] 99%
Re-writing with blocks of size: 256 x 256.
--> [*************************************************************] 99%
Number of valid (not no-data) pixels written: 780916026.
# Finally, generate overviews to enable fast reads at lower resolution
!gdaladdo -r gauss ML_pre_event_3DEP_stereo_bounds.tif 2 4 8 16 32 64 1280...10..]9;4;1;200...10...20 - estimated remaining time: 00:00:1]9;4;1;220...10...20. - estimated remaining time: 00:00:1]9;4;1;250...10...20.. - estimated remaining time: 00:00:1]9;4;1;270...10...20... - estimated remaining time: 00:00:1]9;4;1;300...10...20...30 - estimated remaining time: 00:00:1]9;4;1;320...10...20...30. - estimated remaining time: 00:00:1]9;4;1;350...10...20...30.. - estimated remaining time: 00:00:1]9;4;1;370...10...20...30... - estimated remaining time: 00:00:1]9;4;1;400...10...20...30...40 - estimated remaining time: 00:00:1]9;4;1;420...10...20...30...40. - estimated remaining time: 00:00:1]9;4;1;450...10...20...30...40.. - estimated remaining time: 00:00:1]9;4;1;470...10...20...30...40... - estimated remaining time: 00:00:1]9;4;1;500...10...20...30...40...50 - estimated remaining time: 00:00:1]9;4;1;520...10...20...30...40...50. - estimated remaining time: 00:00:1]9;4;1;550...10...20...30...40...50.. - estimated remaining time: 00:00:1]9;4;1;570...10...20...30...40...50... - estimated remaining time: 00:00:0]9;4;1;600...10...20...30...40...50...60 - estimated remaining time: 00:00:0]9;4;1;620...10...20...30...40...50...60. - estimated remaining time: 00:00:0]9;4;1;650...10...20...30...40...50...60.. - estimated remaining time: 00:00:0]9;4;1;670...10...20...30...40...50...60... - estimated remaining time: 00:00:0]9;4;1;700...10...20...30...40...50...60...70 - estimated remaining time: 00:00:0]9;4;1;720...10...20...30...40...50...60...70. - estimated remaining time: 00:00:0]9;4;1;750...10...20...30...40...50...60...70.. - estimated remaining time: 00:00:0]9;4;1;770...10...20...30...40...50...60...70... - estimated remaining time: 00:00:0]9;4;1;800...10...20...30...40...50...60...70...80 - estimated remaining time: 00:00:0]9;4;1;820...10...20...30...40...50...60...70...80. - estimated remaining time: 00:00:0]9;4;1;850...10...20...30...40...50...60...70...80.. - estimated remaining time: 00:00:0]9;4;1;870...10...20...30...40...50...60...70...80... - estimated remaining time: 00:00:0]9;4;1;900...10...20...30...40...50...60...70...80...90 - estimated remaining time: 00:00:0]9;4;1;920...10...20...30...40...50...60...70...80...90. - estimated remaining time: 00:00:0]9;4;1;950...10...20...30...40...50...60...70...80...90.. - estimated remaining time: 00:00:0]9;4;1;970...10...20...30...40...50...60...70...80...90... - estimated remaining time: 00:00:0]9;4;0;1000...10...20...30...40...50...60...70...80...90...100 - done in 00:00:23.
- Flynn, I. T. W., Bhushan, S., Henderson, S., Corradino, C., Ramsey, M. S., & Collins, E. R. (2026). Satellite data synergy for volcano monitoring: The 2022 Mauna Loa eruption. Journal of Volcanology and Geothermal Research, 474, 108603. 10.1016/j.jvolgeores.2026.108603