globallandcover layer is 'missing'

carlos-mg89 commented 3 years ago

Right now, and as far as I've tested, the overall results are astounding compared to the original OpenMapTiles generator. I'm really happy with the output MBTiles.

However, there's an important layer missing, globallandcover.

The thing is that it isn't part of the OpenMapTiles generator: https://github.com/openmaptiles/openmaptiles

However, if you want to have a great user experience, it's one important layer since it displays grass, forest, tree, scrub and ice for lower zooms (between 0 and 9). Apparently, these aren't included in the output MBTiles, neither in the one from tilemaker, nor in the one from OpenMapTiles project.

So I believe it's because the people from OpenMapTiles decided it was something they didn't want to include. Which is totally fair from their side. They are a company and they have already provided with lots of awesome tools.

I believe it should be possible to generate this layer with the data mentioned above. This is the part related to this layer, in a JSON style:

{
      "id": "globallandcover_grass",
      "type": "fill",
      "paint": {
        "fill-color": "rgba(222, 226, 191, 1)",
        "fill-opacity": {
          "stops": [
            [
              8,
              0.6
            ],
            [
              9,
              0.2
            ]
          ]
        }
      },
      "filter": [
        "all",
        [
          "==",
          "class",
          "grass"
        ]
      ],
      "layout": {
        "visibility": "visible"
      },
      "source": "openmaptiles",
      "maxzoom": 9,
      "minzoom": 0,
      "source-layer": "globallandcover"
    },
    {
      "id": "globallandcover_scrub",
      "type": "fill",
      "paint": {
        "fill-color": "rgba(202, 214, 166, 1)",
        "fill-opacity": 1
      },
      "filter": [
        "all",
        [
          "==",
          "class",
          "scrub"
        ]
      ],
      "layout": {
        "visibility": "visible"
      },
      "source": "openmaptiles",
      "maxzoom": 9,
      "minzoom": 0,
      "source-layer": "globallandcover"
    },
    {
      "id": "globallandcover_tree",
      "type": "fill",
      "paint": {
        "fill-color": "rgba(191, 202, 155, 1)",
        "fill-opacity": 1
      },
      "filter": [
        "all",
        [
          "==",
          "class",
          "tree"
        ]
      ],
      "layout": {
        "visibility": "visible"
      },
      "source": "openmaptiles",
      "maxzoom": 9,
      "minzoom": 0,
      "source-layer": "globallandcover"
    },
    {
      "id": "globallandcover_forest",
      "type": "fill",
      "paint": {
        "fill-color": "rgba(191, 202, 155, 1)",
        "fill-opacity": 1
      },
      "filter": [
        "all",
        [
          "==",
          "class",
          "forest"
        ]
      ],
      "layout": {
        "visibility": "visible"
      },
      "source": "openmaptiles",
      "maxzoom": 9,
      "minzoom": 0,
      "source-layer": "globallandcover"
    },
    {
      "id": "globallandcover_ice",
      "type": "fill",
      "paint": {
        "fill-color": "rgba(255, 255, 255, 1)",
        "fill-opacity": 1
      },
      "filter": [
        "all",
        [
          "==",
          "class",
          "snow"
        ]
      ],
      "layout": {
        "visibility": "visible"
      },
      "source": "openmaptiles",
      "maxzoom": 9,
      "minzoom": 0,
      "source-layer": "globallandcover"
    }

You can observe that the source is the same one, "openmaptiles", but not the layer. The sources that are in the JSON style that MapTiler provides are these ones:

"sources": {
    "contours": {
      "url": "http://www.your-server.local/js/contours.json",
      "type": "vector"
    },
    "landcover": {
      "url": "https://api.maptiler.com/tiles/landcover/tiles.json?key=YOUR_API_KEY",
      "type": "vector"
    },
    "hillshades": {
      "url": "http://www.your-server.local/js/hillshades.json",
      "type": "raster",
      "tileSize": 256
    },
    "openmaptiles": {
      "url": "http://www.your-server.local/js/openmaptiles_base_tiles.json",
      "type": "vector"
    }
  }

The one that includes the PBFs that have the globallandcover layer, are in landcover. The landcover layer doesn't seem to have that data as I've investigated, but I'm not entirely sure about it.

More info here: https://cloud.maptiler.com/tiles/landcover/

P.S. I'll keep researching about this.

systemed commented 3 years ago

Interesting question! Small-scale landcover is often done from external datasets (e.g. CORINE, ESA Landcover, Global Land Cover 2000), which in some cases need to be vectorised from a raster original. It wouldn't be a bad idea to support such a dataset in the same way that we support coastlines.

carlos-mg89 commented 3 years ago

I'm attaching 2 screenshots that represent the following scenarios:

Screenshot with green areas: it uses the maptiler.com 'landcover' source (which has the 'globallandcover' layer)
Screenshot without green areas: it simply doesn't load the maptiler.com source

It's really helpful the information you describe about this layer being generated from external datasets. I'm going to try to find them and see if it's doable from my side to process them.

I find it odd that this isn't possible to be generate it using the OSM data. Since in my opinion, most of the green areas, I believe that they are natural parks, national parks, forests, etc. So I think that it should be possible to achieve a similar result by making these polygons visible in lower zoom levels (between 0 and 9 is where the situation happens).

openmaptiles_with_globallandcover openmaptiles_without_globallandcover

systemed commented 3 years ago

That's a helpful screenshot - thanks!

It should certainly be possible to generate it from OSM data - it's just perhaps not so efficient. Forests are often mapped very 'bittily' in OSM with lots of separate, intricate polygons.

tilemaker can run a simplification algorithm over a single (multi)polygon, but really in a case like this, what you'd want to do is combine the polygons first and then simplify. That's quite CPU-intensive and algorithmically complex. If we can pre-prepare some shapefiles, and then drop them into the config just as we do for the coastline, that would be much simpler.

carlos-mg89 commented 3 years ago

Attached you'd find 2 more screenshots. Both scenarios are using the SAME style.json, with Leaflet. They have the same sources too, except for the base tiles (the ones tilemaker generate).

By the way, I'm using in both scenarios the landcover layer by MapTiler.com, that's why the landcover polygons (natural parks, rivers, glaciers...) are being displayed in both of them.

The information that isn't available in the MBTiles of tilemaker, is the administrative divisions (countries, and states aren't available either in z7-14). Aside from that, some natural areas (in gray in the openmaptiles scenario) aren't present in tilemaker either. Rivers aren't displayed in tilemaker tiles either at this zoom level.

This one, at z5, is using an MBTiles server that I've set up and uses tiles made by tilemaker:

tilemaker_z5

This one, at z5 too, uses an MBTiles generate by OpenMapTiles (https://github.com/openmaptiles/openmaptiles):

openmaptiles_z5

I'm going to start making tests so I can try to make tilemaker to generate tiles that are closer to the OpenMapTiles one. I'll start with the states and countries administrative divisions.

By the way, in zooms 0-4, the ocean disappears in tilemaker tiles, while it remains in openmaptiles tiles.

Tilemaker: tilemaker_z4

OpenMapTiles: openmaptiles_z4

Any help / guidance would be more than welcome.

kleunen commented 3 years ago

Having the administrative divisions would be very nice. Would be great if these can be added.

i have some ideas on combining smaller polygons into larger polygons. Using rtree nearest neigbour can be queried. I will write down some ideas later.

kleunen commented 3 years ago

What you can do is the following, what @systemed suggests for combining the polygons into global polygons is not as difficult as you may think. What you could do is take all the forest polygons and store them in a rtree. Due to the fact that the coordinates of the corners of the polygons are converted to pixels when rendering, precision is lost. Basicly, the coordinates of the polygons are multiplied with some scaling factor and then rounded to a pixel. This causes small polygons to become a single pixel or subpixel at certain zoom level.

This is why the min area and meters -> pixels formula is there: https://github.com/systemed/tilemaker/blob/master/resources/process-openmaptiles.lua#L19-L28 https://github.com/systemed/tilemaker/blob/master/resources/process-openmaptiles.lua#L464-L475

But also, this scaling causes certain polygons to overlap. Small polygons will basicly map to the same pixel. If enough pixels are covered, then actually this pixel should be considered "filled".

The approach I considered is now the following:

Convert the coordinates of all forest areas to pixels at a certain zoom level
Store all the polygons in an rtree
Combine the polygons that overlap. You can go through all polygons and request the nearest polygon. If they overlap, they should be combined.
Combine all the polygons that overlap
Remove all the polygons that do not fit the min.area test
Convert coordinates back to the unzoomed level

This gives the renderer an additional set of polygons which can be used at low zoom levels. I think these polygons should only be calculated for one zoom level. And you only need to run the combination algorithm on polygons which are not visible anymore at the specified zoom level.

kleunen commented 3 years ago

Using OpenStreetMap to Create Land Use and Land Cover Maps:Development of an Application http://pure.iiasa.ac.at/id/eprint/14464/1/See%20chap%207_campelo%202017%20book.pdf

https://github.com/jasp382/glass https://github.com/jasp382/glass/tree/bd820c49c86d05d7cb57afc1a306a4ff988d9c55/core/glass/ete/osm2lulc

It seems the repository is active, so, maybe this person can give some insights

systemed commented 3 years ago

I've added some basic landcover in https://github.com/systemed/tilemaker/commit/40834a06a535f8b44eb50d506f76a07fa0d3152f.

You'll need to create a landcover directory in the same place you're generating your tiles from (typically your tilemaker directory, just as you would for a coastlines directory), and then put in it three sets of unzipped shapefiles from Natural Earth: 10m ice shelves, 10m glaciers, 10m urban areas.

So your directory structure will look like this:

tilemaker
tilemaker/landcover
tilemaker/landcover/ne_10m_urban_areas
tilemaker/landcover/ne_10m_urban_areas/ne_10m_urban_areas.cpg
tilemaker/landcover/ne_10m_urban_areas/ne_10m_urban_areas.dbf

etc.

It'll then use urban areas, glaciers and (if you're generating the Antarctic!) ice shelves in the finished map.

That's all that Natural Earth has, so we'll need to look elsewhere for grassland etc. I'll look at generating some shapefiles from a raster source for this - possibly Global Land Cover, though their download server is currently 500ing.

kleunen commented 3 years ago

But you do not necessairly have to include that in the mbtiles right ? If you use a raster source, you can also combine them in the frontend. If I look at this style: https://github.com/maputnik/osm-liberty

It combines both vector sources as well as raster sources for relief shading. This can work quite well also ?

systemed commented 3 years ago

Yep - for relief shading it's very often the most sensible option (though some people use vectors for that). MBGL has some interesting functionality for actually reading a DEM directly, though I think they missed a trick on compression methods.

For low-zoom landuse we could potentially use the Natural Earth rasters at https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/ .

Shapefiles have more styling flexibility though, so I think it's worth at least trying to generate some. Essentially this would be using gdal_sieve to remove shapes with just a few pixels, then gdal_polygonize to create vectors. I do a few things like this for some of cycle.travel's processing so should be able to have a stab at it.

kleunen commented 3 years ago


REGION=europe

default: all

.PRECIOUS: %.osm.pbf

%.osm.pbf:
    rm -rf $@
    wget --continue http://download.geofabrik.de/$(REGION)/$@ -O $@

landcover/ne_10m_urban_areas/ne_10m_urban_areas.shp:
    mkdir -p landcover/ne_10m_urban_areas
    cd landcover/ne_10m_urban_areas && rm -rf ne_10m_urban_areas.zip
    cd landcover/ne_10m_urban_areas && wget https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_urban_areas.zip
    cd landcover/ne_10m_urban_areas && unzip -o ne_10m_urban_areas.zip

landcover/ne_10m_antarctic_ice_shelves_polys/ne_10m_antarctic_ice_shelves_polys.shp:
    mkdir -p landcover/ne_10m_antarctic_ice_shelves_polys
    cd landcover/ne_10m_antarctic_ice_shelves_polys && rm -rf ne_10m_antarctic_ice_shelves_polys.zip
    cd landcover/ne_10m_antarctic_ice_shelves_polys && wget https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_antarctic_ice_shelves_polys.zip
    cd landcover/ne_10m_antarctic_ice_shelves_polys && unzip -o ne_10m_antarctic_ice_shelves_polys.zip

landcover/ne_10m_glaciated_areas/ne_10m_glaciated_areas.shp:
    mkdir -p landcover/ne_10m_glaciated_areas
    cd landcover/ne_10m_glaciated_areas rm -rf ne_10m_glaciated_areas.zip
    cd landcover/ne_10m_glaciated_areas && wget https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_glaciated_areas.zip
    cd landcover/ne_10m_glaciated_areas && unzip -o ne_10m_glaciated_areas.zip

coastline/water_polygons.shp:
    mkdir -p coastline
    cd coastline && rm -rf water-polygons-split-4326.zip
    cd coastline && wget https://osmdata.openstreetmap.de/download/water-polygons-split-4326.zip
    cd coastline && unzip -o water-polygons-split-4326.zip
    cd coastline && mv water-polygons-split-4326/* .

%.mbtiles: %.osm.pbf coastline/water_polygons.shp landcover/ne_10m_glaciated_areas/ne_10m_glaciated_areas.shp landcover/ne_10m_antarctic_ice_shelves_polys/ne_10m_antarctic_ice_shelves_polys.shp landcover/ne_10m_urban_areas/ne_10m_urban_areas.shp
    ./tilemaker --init-store 500:100 --input $< --output $@ --config config-openmaptiles.json --process process-openmaptiles.lua

all: liechtenstein-latest.mbtiles

systemed / tilemaker

globallandcover layer is 'missing' #203