gsi-dem-toolkit 0.1.0

Download GSI terrain tiles, decode proprietary PNG elevation format, and compute slope/TRI/Tobler hiking cost for Japan

gsi-dem-toolkit

Download GSI terrain tiles, decode Japan's proprietary PNG elevation format, and compute slope/TRI/Tobler hiking cost.

Install

pip install gsi-dem-toolkit

Quick Start

from gsi_dem_toolkit import run_pipeline, create_visualization

# Download DEM, compute terrain metrics, output GeoParquet
output = run_pipeline(prefecture="tokyo", zoom=14)

# Generate interactive HTML map
create_visualization(output)
print(f"Results: {output}")
print(f"Map: {output.with_suffix('.html')}")

What is GSI DEM?

The Geospatial Information Authority of Japan (GSI) provides free 5m and 10m resolution Digital Elevation Model tiles covering all of Japan. The elevation data is encoded in a proprietary RGB PNG format where each pixel stores centimeter-precision elevation as a 24-bit integer across the R, G, and B channels. This toolkit handles tile download, PNG decoding, GeoTIFF mosaicking, and terrain analysis -- turning raw GSI tiles into actionable slope, ruggedness, and walkability metrics.

API Reference

Pipeline

run_pipeline(prefecture="tokyo", bbox=None, zoom=14, layer="dem5a_png",
             cell_size_deg=0.005, zones_file=None, output_dir=None,
             max_workers=4) -> Path

Full pipeline: download DEM tiles, compute slope/TRI/Tobler rasters, create analysis grid, compute zonal statistics, calculate Terrain Difficulty Score. Returns path to output GeoParquet.

Tile Operations

download_dem_bbox(bbox, zoom=14, layer="dem5a_png", output_path=None,
                  max_workers=4) -> tuple[ndarray, dict]

Download and mosaic all DEM tiles for a bounding box into a single raster. Optionally saves as GeoTIFF. Returns (elevation_array, rasterio_profile).

decode_dem_png(png_bytes) -> ndarray

Decode a single GSI PNG tile to a 256x256 float32 elevation array in meters.

lat_lon_to_tile(lat, lon, zoom) -> tuple[int, int]

Convert WGS84 coordinates to tile (x, y) indices at a given zoom level.

tile_bounds(x, y, zoom) -> tuple[float, float, float, float]

Get the WGS84 bounding box (west, south, east, north) for a tile.

get_tiles_for_bbox(bbox, zoom) -> list[tuple[int, int]]

Get all tile (x, y) indices covering a bounding box.

Terrain Analysis

compute_slope(dem, transform, crs="EPSG:4326") -> ndarray

Compute slope in degrees using Horn's (1981) 3x3 gradient method.

compute_tri(dem) -> ndarray

Compute Terrain Ruggedness Index using Riley et al. (1999) method.

compute_tobler_cost(slope_deg) -> ndarray

Compute Tobler hiking function cost in minutes per kilometer.

Aggregation

compute_zonal_stats(zones, slope_raster, tri_raster) -> GeoDataFrame

Aggregate raster metrics to polygon zones. Adds columns: slope_mean, slope_max, slope_std, slope_median, slope_p90, pct_steep_5deg, pct_steep_10deg, tobler_cost_min_per_km, tri_mean, tri_max.

compute_terrain_difficulty_score(gdf, weights=None) -> GeoDataFrame

Composite Terrain Difficulty Score (0=flat, 1=extreme). Adds terrain_difficulty_score and terrain_class (Flat/Mild/Moderate/Difficult/Very Difficult) columns.

create_grid_zones(bbox, cell_size_deg=0.005) -> GeoDataFrame

Generate a regular grid of square polygons for analysis. Default ~500m cells.

Visualization

create_visualization(parquet_path, output_html=None) -> Path

Create an interactive Folium choropleth map of Terrain Difficulty Score with GSI base map layers.

Reference Data

PREFECTURE_BBOX: dict[str, tuple[float, float, float, float]]

WGS84 bounding boxes (west, south, east, north) for all 47 Japanese prefectures.

GSI PNG Encoding

GSI encodes elevation in the RGB channels of a 256x256 PNG tile:

raw_value = R * 65536 + G * 256 + B
if raw_value < 2^23:    elevation_m = raw_value * 0.01
if raw_value == 2^23:   no data (NaN)
if raw_value > 2^23:    elevation_m = (raw_value - 2^24) * 0.01   # negative elevation

Resolution depends on the DEM layer: dem5a_png (5m mesh A, best), dem5b_png, dem5c_png, dem_png (10m, full coverage fallback). The toolkit automatically falls back through layers if the preferred one lacks data for a tile.

Data Sources

• Elevation tiles: GSI (cyberjapandata.gsi.go.jp)
• DEM layers: dem5a_png (5m, best resolution) through dem_png (10m, full coverage)
• Base maps: GSI Pale and Relief tile layers (used in visualization)

License

MIT -- see LICENSE.