Finland Forest 2019 Creation

From OpenStreetMap Wiki
(Redirected from Posiki python script)
Jump to navigation Jump to search

This page describe how Finland Forest coverage import was attempted.

Background information

Finnish Forest cover has been problematic since 2012 (or even earlier). Check more details from Finland Corine 2006 Import

Process description

Brief process workflow

  1. Forest are created from LUKE's multi-source national forest inventory (MS-NFI) from 2015 (download service, license)
  2. Clip out OSM features from forest plots
  3. Clip out National Land Survey of Finland's Topographic Database features from forest plots
  4. Import forest plots to OSM

Data sources

All raw source data sets are openly licensed.

Process workflow

Note: all processes are based on Finnish national map sheet grid. This is most convient way to handle Finnish GIS datasets, you download grid from NLS.FI's download service.

Load MS_NFI data to PostGIS

Loading of MS_NFI raster data to PostGIS is done with raster2pgsql -command. Basic command is:

raster2pgsql -s 3067 -t 100x100 -a tilavuus_vmi1x_1216_Q4.tif mvmi.tilavuus

Finnish coordinate system is EPSG:3067, we use "volume of growing stock" datasaset (Puuston tilavuus kuutioina (in Finnish)). LUKE is delivering data in 200k map sheets.

MS_NFI raster data: green forest volume is more than 3 cubicmeters, pink less than 3 and white: no forest.

Raw source data has information about volume of forest in cubic meters, cell size is 16 x 16 meters.

MS_NFI raw forestry information (16 x 16 m pixel size) with NLS.fi ortophoto

Vectorize forestry data

Next step is to vectorize forestry data with PostGIS ST_DumpAsPolygons -function. Following SQL-pseudo-code also includes selection of data for 25k map sheets (for easier/faster processing):

CREATE TABLE q4233.forest_q4233 AS (
	WITH aoi AS (
	 SELECT geom, lehtitunnu
	 FROM nlsfi.utm25
	 WHERE lehtitunnu = upper('q4233')
	)

	SELECT val,geom
	FROM (
		SELECT (ST_DumpAsPolygons(a.rast)).*  
		FROM mvmi.tilavuus a, aoi b
		WHERE ST_Intersects(a.rast, b.geom)
		) AS foo
	WHERE val > 3
);
Vectorised forest squares

Avoid sharp lines

To avoid sharp lines, we create circles around squares, with PostGIS function ST_MinimumBoundingCircle-function:

UPDATE q4233.forest_q4233 
	SET geom_circle = ST_MinimumBoundingCircle(geom);
Minimum bounding circles around forest squares

Clustering forests

Then we cluster forestry circles with ST_ClusterIntersecting -function:

...
clust AS (
 SELECT unnest(ST_ClusterIntersecting(geom_circle)) as geom
 FROM forest
...

Note: ST_ClusterIntersecting has limitations (PostgreSQL arrays size is limiting), so all process after this is made with 10k sheets (6 km x 6 km squares).

Clustered forests, also smaller areas (less than 1000 square meters are removed)

Smoothing borders

Next step will be smoothing borders with ST_ChaikinSmoothing -function:

CREATE TABLE q4233.forest_chs_q4233d AS (
  SELECT ST_MakeValid(ST_ChaikinSmoothing(ST_SimplifyPreserveTopology(geom,12),5)) as geom
  FROM q4233.forest_plot_q4233d
);
Smoothed forests.png

Comparing to existing OpenStreetMap features

This creation will not overwrite any existing natural or landuse tagged features. This migh be very conservative way to handle automatic mapping, but this is how I decide to handle situation.

Also highway -features are buffered with following buffered sizes (meters):

... 
CASE
  WHEN highway = 'motorway' THEN 10
  WHEN highway = 'trunk' THEN 8
  WHEN highway = 'primary' THEN 8
  WHEN highway = 'secondary' THEN 7
  WHEN highway = 'tertiary' THEN 5 
  WHEN highway = 'unclassified' THEN 4 
  WHEN highway = 'service' THEN 4 
  WHEN highway = 'residential' THEN 4 
  WHEN highway = 'track' THEN 3 
  WHEN highway = 'cycleway' THEN 2 
  WHEN highway = 'footway' THEN 2 
  WHEN highway = 'path' THEN 1 
  ELSE 0.5 
END
...
Clip With OSM features.png

Comparing to National Land Survey of Finland's features

For preparing future and ongoing imports of NLS.FI's data sets, forestry plots are clipped also with NLS.FI topographic database features.

Area features
NLS.fi name in English Notes
Maatalousmaa Farmland
Niitty Meadow
Virtavesialue River (area)
Järvi Lake
Line features
NLS.fi name in English Notes
Sähkölinja Electric lines Class 22311 (high voltage) buffered with 20 meters

Class 22312 (low voltage) buffered with 5 meters

Tieviiva Highways Class 12111 buffered with 10 meters

Class 12112 buffered with 8 meters

Class 12121 buffered with 6 meters

Class 12122 buffered with 6 meters

Class 12131 buffered with 3 meters

Class 12132 buffered with 3 meters

Class 12141 buffered with 2 meters

Class 12313 buffered with 0,5 meters

Class 12314 buffered with 1,5 meters

Class 12316 buffered with 1,5 meters


Green polygons are clipped forest polygons, other colors are NLS.FI's features.

Finalisation of forest polygons

Forestry polygons are finally slighly simplifies with ST_SimplifyPreserveTopology:

UPDATE q4233.forest_final_q4233d 
	SET geom = ST_Force2D(ST_SimplifyPreserveTopology(geom,1));

Also smallish polygons (less than 1000 m²) are deleted.

OSM Final Forest.gif
Major phases in creation of forestry polygons

Preparation for loading to OpenStreetMap

Limit of OSM changeset is 10k elements. Separated PostGIS function has been developed to group forestry polygons to at maximum 9000 points / changeset.

ogr2osm.py has been used to create OSM import files.

Loading to OpenStreetMap

Modified version of bulk_upload.py has been used. Modification was mainly convert to use Python 3 ( orginal version of bulk_upload.py is Python 2 script).

Loading has been done with posiki_import_forest -account. Changeset has been named with 10K sheet names.

Sample area has been uploaded to OSM with two (2) changesets:

Problems

Self-intersecting forestry plots

Easiest way to detect self-intersecting forestry polygons is Geofabrik's OSM Inspector.

Self-intersecting polygon.png

Solution

JOSM is easy tool to correct these kind of problems. Just download polygon to JOSM and delete one point. Sometimes self-intersections could be triggier to correct.

JOSM with self-intersecting polygon

Disclaimer: in PostGIS this is valid polygon, but not in OpenStreetMap. This is how GeoWorld is different in differenent spaces ;-)


Removing forest plots

After feedback and discussion in OSM Forum, I have decided to remove all forestry imports which are done by me.

Method

All created forest plots (nodes, ways and relations) will be removed, unless other users has _not_ edited features.

Progress

Reverting of forests has been started on May 2021 and end November 2021

Misc notes

Feedback

Send by email to pos (at) iki.fi or via by OSM account (posiki_import_forest or posiki)

Disclaimer

This is not really import, but automatic creation of new features from open data sources. IMHO. I agree your disagree, if you disagree.