AW3D30 elevation over Mt. Fuji — research scenario¶

Pull ALOS World 3D 30 m DSM for the full Mt. Fuji caldera + foothills (0.5° × 0.5° bbox, ~55 × 55 km) through the authless jaxa-earth protocol, then derive elevation stats and a representative profile. This is the kind of pull a glaciologist or volcanologist runs to bound a study area before higher-resolution fieldwork.

The fetch uses the official jaxa.earth API; the COG write goes through pyramids. Needs the [jaxa] extra; no credentials.

Setup¶

In [ ]:

import tempfile
from pathlib import Path

import numpy as np

from earthlens.core import EarthLens

import tempfile from pathlib import Path import numpy as np from earthlens.core import EarthLens

Request¶

Mt. Fuji summit sits at ~35.36° N, 138.73° E (peak elevation 3776 m). Bbox the caldera + foothills (lat 35.1 – 35.6, lon 138.5 – 139.0) and pull at 1000 px/° (~110 m pixels). AW3D30 is a static product — the SDK snaps to its single available epoch (2021-02) once a wide date window is passed.

In [ ]:

out_dir = Path(tempfile.mkdtemp(prefix='earthlens-jaxa-aw3d30-'))
lens = EarthLens(
    data_source='jaxa',
    variables=['elevation'],
    start='2000-01-01',
    end='2030-12-31',
    lat_lim=[35.10, 35.60],
    lon_lim=[138.50, 139.00],
    resolution=1000.0,
    path=out_dir,
)

out_dir = Path(tempfile.mkdtemp(prefix='earthlens-jaxa-aw3d30-')) lens = EarthLens( data_source='jaxa', variables=['elevation'], start='2000-01-01', end='2030-12-31', lat_lim=[35.10, 35.60], lon_lim=[138.50, 139.00], resolution=1000.0, path=out_dir, )

Download¶

In [ ]:

written = lens.download()
for p in written:
    print(f'{p.name}  ({p.stat().st_size // 1024} KB)')

written = lens.download() for p in written: print(f'{p.name} ({p.stat().st_size // 1024} KB)')

Open the COG + derive elevation statistics¶

Pyramids reads it back in EPSG:4326. We compute the summary statistics you'd quote in a study-area description: min/max/mean elevation, the peak's geographic location, and a rough hypsometric breakdown.

In [ ]:

from pyramids.dataset import Dataset

cog = Dataset.read_file(str(written[0]))
geo = cog.geotransform
arr = cog.read_array().astype('float32')

print(f'EPSG:        {cog.epsg}')
print(f'array shape: {arr.shape}  ({arr.size:,} pixels)')
print(f'min elev:    {arr.min():.0f} m')
print(f'max elev:    {arr.max():.0f} m  (Mt. Fuji peak is 3776 m)')
print(f'mean elev:   {arr.mean():.0f} m')
print(f'pct > 2500m: {(arr > 2500).mean() * 100:.1f} %  (above tree line)')
print(f'pct > 3500m: {(arr > 3500).mean() * 100:.1f} %  (summit cone)')

from pyramids.dataset import Dataset cog = Dataset.read_file(str(written[0])) geo = cog.geotransform arr = cog.read_array().astype('float32') print(f'EPSG: {cog.epsg}') print(f'array shape: {arr.shape} ({arr.size:,} pixels)') print(f'min elev: {arr.min():.0f} m') print(f'max elev: {arr.max():.0f} m (Mt. Fuji peak is 3776 m)') print(f'mean elev: {arr.mean():.0f} m') print(f'pct > 2500m: {(arr > 2500).mean() * 100:.1f} % (above tree line)') print(f'pct > 3500m: {(arr > 3500).mean() * 100:.1f} % (summit cone)')

Locate the peak in geographic coordinates¶

Convert the argmax pixel back to (lat, lon) using the geotransform's north-up convention: (lon0, dx, 0, lat0, 0, -dy).

In [ ]:

row, col = np.unravel_index(arr.argmax(), arr.shape)
lon0, dx, _, lat0, _, neg_dy = geo
lat = lat0 + neg_dy * row  # neg_dy is negative
lon = lon0 + dx * col
print(f'peak pixel: row={row}, col={col}, elevation={arr.max():.0f} m')
print(f'geographic: lat={lat:.4f} N, lon={lon:.4f} E')
print('Mt. Fuji official: lat=35.3606 N, lon=138.7274 E')

row, col = np.unravel_index(arr.argmax(), arr.shape) lon0, dx, _, lat0, _, neg_dy = geo lat = lat0 + neg_dy * row # neg_dy is negative lon = lon0 + dx * col print(f'peak pixel: row={row}, col={col}, elevation={arr.max():.0f} m') print(f'geographic: lat={lat:.4f} N, lon={lon:.4f} E') print('Mt. Fuji official: lat=35.3606 N, lon=138.7274 E')

Hypsometric histogram¶

Distribution of elevations across the bbox — useful for picking a reasonable contour interval or planning DEM-conditioned hydrology.

In [ ]:

edges = list(range(0, 4001, 250))
counts, _ = np.histogram(arr.ravel(), bins=edges)
for lo, hi, n in zip(edges[:-1], edges[1:], counts):
    bar = '#' * (n * 40 // max(counts.max(), 1))
    print(f'{lo:>4}-{hi:<4} m: {n:>6}  {bar}')

edges = list(range(0, 4001, 250)) counts, _ = np.histogram(arr.ravel(), bins=edges) for lo, hi, n in zip(edges[:-1], edges[1:], counts): bar = '#' * (n * 40 // max(counts.max(), 1)) print(f'{lo:>4}-{hi:<4} m: {n:>6} {bar}')