# 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
```

## Public API

```python
from gsi_dem_toolkit import (
    run_pipeline,
    download_dem_bbox,
    decode_dem_png,
    lat_lon_to_tile,
    tile_bounds,
    get_tiles_for_bbox,
    compute_slope,
    compute_tri,
    compute_tobler_cost,
    compute_zonal_stats,
    compute_terrain_difficulty_score,
    create_grid_zones,
    create_visualization,
    PREFECTURE_BBOX,
)
```

### 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 -> slope/TRI/Tobler rasters -> zonal stats -> Terrain Difficulty Score -> GeoParquet.

### download_dem_bbox(bbox, zoom=14, layer="dem5a_png", output_path=None, max_workers=4) -> tuple[ndarray, dict]
Download + mosaic DEM tiles to GeoTIFF. Returns (elevation_array, rasterio_profile).

### decode_dem_png(png_bytes) -> ndarray
Decode single GSI PNG tile (256x256) to float32 elevation array in meters.

### lat_lon_to_tile(lat, lon, zoom) -> tuple[int, int]
WGS84 coordinates to tile (x, y) indices.

### tile_bounds(x, y, zoom) -> tuple[float, float, float, float]
Tile indices to WGS84 bbox (west, south, east, north).

### get_tiles_for_bbox(bbox, zoom) -> list[tuple[int, int]]
All tile indices covering a bounding box.

### compute_slope(dem, transform, crs="EPSG:4326") -> ndarray
Horn's method slope in degrees. Handles degree-based CRS by converting cell size to meters.

### compute_tri(dem) -> ndarray
Riley et al. Terrain Ruggedness Index in meters.

### compute_tobler_cost(slope_deg) -> ndarray
Tobler hiking function: minutes per kilometer.

### compute_zonal_stats(zones, slope_raster, tri_raster) -> GeoDataFrame
Aggregate raster metrics to polygons. Adds slope_mean, slope_max, slope_std, 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 score 0-1. Adds terrain_difficulty_score and terrain_class columns.

### create_grid_zones(bbox, cell_size_deg=0.005) -> GeoDataFrame
Generate regular grid polygons (~500m cells) for analysis zones.

### create_visualization(parquet_path, output_html=None) -> Path
Folium choropleth map of Terrain Difficulty Score with GSI base maps.

### PREFECTURE_BBOX: dict[str, tuple[float, float, float, float]]
Bounding boxes for all 47 prefectures. Keys are lowercase romaji.

## Common Pattern

```python
from gsi_dem_toolkit import run_pipeline, create_visualization

output = run_pipeline(prefecture="osaka", zoom=14)
create_visualization(output)
```

## Key Details
- GSI PNG encoding: R*65536 + G*256 + B, threshold at 2^23 for no-data, >2^23 for negative elevation
- DEM layers: dem5a_png (5m best) -> dem5b_png -> dem5c_png -> dem_png (10m fallback)
- Cache: ~/.cache/gsi-dem-toolkit/ (raw tiles, processed rasters)
- Rate limit: 50ms delay between tile requests; be respectful to GSI servers
- Requires Python >= 3.10, rasterio, rasterstats, numpy, geopandas, Pillow
