h3-landcover

Make low resolution landcover vector tiles from high resolution raster landcover data with h3 downsampling

Demo

https://wipfli.github.io/h3-landcover/

Design

The following design descriptions for the individual steps were given to ChatGPT which then produced the python scripts. I had to tweek some of the python scripts by hand afterwards, but the majority of the code content was written automatically by ChatGPT.

Download

Download raster landcover data from https://lcviewer.vito.be/download. The data is split into 20 by 20 degree tiles. Use 2019 version, 100 m resolution.

Attribution landcover data: https://lcviewer.vito.be/about

The download URLs look like:

https://s3-eu-west-1.amazonaws.com/vito.landcover.global/v3.0.1/2019/E000N60/E000N60_PROBAV_LC100_global_v3.0.1_2019-nrt_Bare-CoverFraction-layer_EPSG-4326.tif

template:

https://s3-eu-west-1.amazonaws.com/vito.landcover.global/v3.0.1/2019/{longitude}{latitude}/{longitude}{latitude}_PROBAV_LC100_global_v3.0.1_2019-nrt_{kind}-CoverFraction-layer_EPSG-4326.tif

latitudes: S40, S20, N00, N20, N40, N60, N80 longitudes: W180, W160, W140, W120, W100, W080, W060, W040, W020, E000, E020, E040, E060, E080, E100, E120, E140, E160 kinds: Snow, Bare, Grass, Crops, Tree, BuiltUp

Download all files into a sources folder.

Write a python script which loops through the filenames and downloads them to the sources folder with system calls using the cli tool axel which has parallel downloads with the -a -n 10 option.

download.py

The total landcover raster data is 26 GB in 564 files.

Raster to H3 Classification

For every image in the sources folder:

• Open the image
• Loop through the pixels
• The pixel value is landcover in percent, the pixel size is 100 m by 100 m
• Call ReadAsArray only once per image
• Make sure to ignore nodata pixels
• Print fractional progress after every 10000th pixel
• Multiply pixel area by pixel value divided by 100 to the the covered area of the pixel
• Hold a hash map H3 index to covered area
• Add the area to the H3 cell in the hash map in which the pixel is located, use H3 resolution 9
• Save the hashmap as a json file. Save it at classified/{original_filename}.json and create the folder if needed.
• Use the osgeo gdal library
• Use h3.api.numpy_int otherwise it is too slow

Write a python script for this.

Second promt: Using the above script, write a version which runs on a 12 core CPU. The images are independent and can be treated in parallel.

classify.py

Combine Grass, Bare, Snow

The kinds Grass, Bare, Snow are plotted on top of eachother. In mountains, we first have Grass landcover, then with higher altitude, Bare landcover, and above this we have Snow. To not have any gaps between Snow and Bare landcover, lets combine the following layers:

• Grass-Bare-Snow
• Bare-Snow
• Snow

The Snow layer we have already. Bare-Snow, and Grass-Bare-Snow have to be generated.

The folder classified contains files with template filename {longitude}{latitude}_PROBAV_LC100_global_v3.0.1_2019-nrt_{kind}-CoverFraction-layer_EPSG-4326.json.

With

latitudes: S40, S20, N00, N20, N40, N60, N80 longitudes: W180, W160, W140, W120, W100, W080, W060, W040, W020, E000, E020, E040, E060, E080, E100, E120, E140, E160

The file content has the form

{
  "1453143554": 4636.3,
  "1453459834": 634.0
}

i.e. dictionary with key H3 cell id, and value covered area of cell in m2.

Loop through all latitude, longitude combinations. If the file exists, loop through the kinds [Bare, Snow] and [Grass, Bare, Snow]. Create a new dictionnary, loop through the cells of the input dictionary and add them to the cells of the new output dictionarry.

Store the resulting H3 map in a file called classified/{longitude}{latitude}_PROBAV_LC100_global_v3.0.1_2019-nrt_{new_kind}-CoverFraction-layer_EPSG-4326.json where new_kind is Bare-Snow or Grass-Bare-Snow.

combine.py

Size of classified folder: 271 GB

H3 to GeoJSON

The folder classified contains files with json H3 maps. The file contents are of the form

{
  "1453143554": 4636.3,
  "1453459834": 634.0
}

i.e. H3 index as key and covered area in m2 as value. The files are in H3 resolution 9. For every file, do the following:

• Read the original file contents to input_resolution_map variable.
• Loop through the following target resolutions 9, 8, 7, 6, 5, 4
• Downsample if the target resolution is not 9:
  • Create an empty target_resolution_map.
  • Loop through the entries of the input_resolution_map
  • For every input H3 index at resolution 9, find the parent H3 index at the target resolution
  • Add the input_resolution_map value to the target_resolution_map at the parent H3 index
• If the target resolution is 9, use the the input_resolution_map as the target_resolution map
• Thresholding: For every cell in the target_resolution_map, decide if more than 50 percent of the cell area is covered. If a cell is covered, add its H3 id to a list called activated.
• Convert to GeoJSON: Convert the activated H3 cell list to a GeoJSON multipolygon.
• Store the multipolygon at the following template filename geojson/{original_filename}-resolution-{target_resolution}.geojson.

to_geojson.py

Size geojson folder: 49 GB

Arctic and Antarctic Polygons

The Arctic and Antarctic are outside of the coverage region of the landcover data. Create polygons for these regions separately which show that they are Snow.

Steps:

• Create a geojson polygons with the bounds South 86 deg to 59 deg, E 180 deg to W 180 deg
• For the resolutions 9, 8, 7, 6, 5, 4, save this polygon in the following location geojson/ANTARCTIC_PROBAV_LC100_global_v3.0.1_2019-nrt_Snow-CoverFraction-layer_EPSG-4326.json-resolution-{resolution}.geojson
• Create a geojson polygons with the bounds North 77 deg to 86 deg, E 180 deg to W 180 deg
• For the resolutions 9, 8, 7, 6, 5, 4, save this polygon in the following location geojson/ARCTIC_PROBAV_LC100_global_v3.0.1_2019-nrt_Snow-CoverFraction-layer_EPSG-4326.json-resolution-{resolution}.geojson

Write a python script for this without using a GeoJSON library. Just write json directly.

arctic_antarctic.py

GeoPackage Zip Files

To consume the geometries in Planetiler, we need to transform the GeoJSON files to GeoPackage files. Files of the same kind should be bundled in a zip file.

Steps are:

• Create the zips folder if needed
• Loop through the kinds Snow, Bare-Snow, Grass-Bare-Snow, Crops, Tree, BuiltUp
• Loop through the resolutions 9, 8, 7, 6, 5, 4
• For all files in the geojson folder which match the pattern *_PROBAV_LC100_global_v3.0.1_2019-nrt_{kind}-CoverFraction-layer_EPSG-4326.json-resolution-{resolution}.geojson do:
  • Convert the file from GeoJSON to GPKG
  • Add the file to a zip folder
  • Save the Zip folder at zips/{kind}-resolution-{resolution}.zip

Write a python script to do this. You can assume that the ogr2ogr utility is installed.

gpkg_zip.py

Size zips folder: 13 GB

Planetiler

Run planetiler with

cd planetiler
./regenerate-swissmap.sh

Download

https://pub-726b01260c98468a9387cc0dfcb7386b.r2.dev/h3-landcover.pmtiles (282 MB)