h3-landcover

April 21, 2023 · View on GitHub

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

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