Metadata-Version: 2.4
Name: gsi-dem-toolkit
Version: 0.1.0
Summary: Download GSI terrain tiles, decode proprietary PNG elevation format, and compute slope/TRI/Tobler hiking cost for Japan
Project-URL: Homepage, https://github.com/Leuvtern/gsi-dem-toolkit
Project-URL: Repository, https://github.com/Leuvtern/gsi-dem-toolkit
Project-URL: Issues, https://github.com/Leuvtern/gsi-dem-toolkit/issues
Author-email: William Ma <william.ma.japan@gmail.com>
License-Expression: MIT
License-File: LICENSE
Keywords: dem,elevation,geospatial,gis,gsi,japan,slope,terrain
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Developers
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: GIS
Requires-Python: >=3.10
Requires-Dist: branca>=0.7
Requires-Dist: folium>=0.15
Requires-Dist: geopandas>=0.14
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: pillow>=10.0
Requires-Dist: rasterio>=1.3
Requires-Dist: rasterstats>=0.19
Requires-Dist: requests>=2.31
Requires-Dist: shapely>=2.0
Description-Content-Type: text/markdown

# gsi-dem-toolkit

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

## Install

```bash
pip install gsi-dem-toolkit
```

## Quick Start

```python
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](https://www.gsi.go.jp/) (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

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

```python
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)`.

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

```python
lat_lon_to_tile(lat, lon, zoom) -> tuple[int, int]
```
Convert WGS84 coordinates to tile (x, y) indices at a given zoom level.

```python
tile_bounds(x, y, zoom) -> tuple[float, float, float, float]
```
Get the WGS84 bounding box `(west, south, east, north)` for a tile.

```python
get_tiles_for_bbox(bbox, zoom) -> list[tuple[int, int]]
```
Get all tile (x, y) indices covering a bounding box.

### Terrain Analysis

```python
compute_slope(dem, transform, crs="EPSG:4326") -> ndarray
```
Compute slope in degrees using Horn's (1981) 3x3 gradient method.

```python
compute_tri(dem) -> ndarray
```
Compute Terrain Ruggedness Index using Riley et al. (1999) method.

```python
compute_tobler_cost(slope_deg) -> ndarray
```
Compute Tobler hiking function cost in minutes per kilometer.

### Aggregation

```python
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.

```python
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.

```python
create_grid_zones(bbox, cell_size_deg=0.005) -> GeoDataFrame
```
Generate a regular grid of square polygons for analysis. Default ~500m cells.

### Visualization

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

```python
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](https://maps.gsi.go.jp/development/ichiran.html) (`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](LICENSE).
