Merge pull request #294 from statisticsnorway/map-hatches

Wms with polygon masks

Author: mortewle

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "ssb-sgis"
-version = "1.2.7"
+version = "1.2.8"
 description = "GIS functions used at Statistics Norway."
 authors = ["Morten Letnes <morten.letnes@ssb.no>"]
 license = "MIT"

src/sgis/geopandas_tools/buffer_dissolve_explode.py

@@ -432,7 +432,7 @@ def _run_func_by_cluster(
 
     def get_group_clusters(group: GeoDataFrame):
         """Adds cluster column. Applied to each group because much faster."""
-        group = group.reset_index(drop=True)
+        group = group[~group.is_empty].reset_index(drop=True)
         group["_cluster"] = get_cluster_mapper(group, predicate=predicate)
         group["_cluster"] = get_grouped_centroids(group, groupby="_cluster")
         return group

src/sgis/maps/explore.py

@@ -779,9 +779,13 @@ def _add_tiles(
         for tile in tiles:
             to_tile(tile, max_zoom=self.max_zoom).add_to(mapobj)
 
-    def _add_wms(self, map_: folium.Map, bbox: Any) -> None:
+    def _add_wms(self, map_: folium.Map) -> None:
+        try:
+            gdf = self._gdf.to_crs(4326)
+        except Exception:
+            gdf = self._gdf.set_crs(4326)
         for wms in self.wms:
-            tiles = wms.get_tiles(bbox, max_zoom=self.max_zoom)
+            tiles = wms.get_tiles(gdf, max_zoom=self.max_zoom)
             for tile in tiles.values():
                 map_.add_child(tile)
 
@@ -995,7 +999,7 @@ def _make_folium_map(
             # folium.LayerControl(collapsed=False).add_to(m)
 
         if self.wms:
-            self._add_wms(m, bounds)
+            self._add_wms(m)
 
         return m
 

src/sgis/maps/maps.py

@@ -197,10 +197,13 @@ def explore(
 
         bounds: Polygon = box(*get_total_bounds(*gdfs, *list(kwargs.values())))
 
-        any_intersections: bool = mask.intersects(bounds).any()
+        if bounds is not None:
+            any_intersections: bool = mask.intersects(bounds).any()
+        else:
+            any_intersections = False
         if not any_intersections and to_crs is None:
             mask = to_gdf(Polygon(), to_crs)
-        elif not any_intersections:
+        elif not any_intersections and bounds is not None:
             bounds4326 = to_gdf(bounds, to_crs).to_crs(25833).geometry.iloc[0]
             mask4326 = mask.set_crs(4326, allow_override=True).to_crs(25833)
 

src/sgis/maps/norge_i_bilder.json

@@ -4,18 +4,33 @@
 import re
 from collections.abc import Iterable
 from dataclasses import dataclass
+from io import BytesIO
 from pathlib import Path
 from typing import Any
 from urllib.request import urlopen
 
 import folium
+import numpy as np
+import pandas as pd
 import shapely
-
+from geopandas import GeoDataFrame
+from geopandas import GeoSeries
+from shapely import Geometry
+from shapely import get_exterior_ring
+from shapely import make_valid
+from shapely import polygons
+from shapely import set_precision
+from shapely import simplify
+from shapely.errors import GEOSException
+
+from ..geopandas_tools.conversion import to_gdf
 from ..geopandas_tools.conversion import to_shapely
+from ..geopandas_tools.sfilter import sfilter
+from ..raster.image_collection import Band
 
 JSON_PATH = Path(__file__).parent / "norge_i_bilder.json"
 
-JSON_YEARS = [str(year) for year in range(1999, 2025)]
+JSON_YEARS = [str(year) for year in range(2006, datetime.datetime.now().year + 1)]
 
 DEFAULT_YEARS: tuple[str] = tuple(
     str(year)
@@ -57,7 +72,7 @@ class NorgeIBilderWms(WmsLoader):
     show: bool | Iterable[int] | int = False
     _use_json: bool = True
 
-    def load_tiles(self) -> None:
+    def load_tiles(self, verbose: bool = False) -> None:
         """Load all Norge i bilder tiles into self.tiles."""
         url = "https://wms.geonorge.no/skwms1/wms.nib-prosjekter?SERVICE=WMS&REQUEST=GetCapabilities"
 
@@ -107,7 +122,10 @@ def load_tiles(self) -> None:
                 this_tile["name"] = name
                 this_tile["bbox"] = this_bbox
                 year = name.split(" ")[-1]
-                if year.isnumeric() and len(year) == 4:
+                is_year_or_interval: bool = all(
+                    part.isnumeric() and len(part) == 4 for part in year.split("-")
+                )
+                if is_year_or_interval:
                     this_tile["year"] = year
                 else:
                     this_tile["year"] = "9999"
@@ -116,41 +134,123 @@ def load_tiles(self) -> None:
 
         self.tiles = sorted(all_tiles, key=lambda x: x["year"])
 
-    def get_tiles(self, bbox: Any, max_zoom: int = 40) -> list[folium.WmsTileLayer]:
-        """Get all Norge i bilder tiles intersecting with a bbox."""
+        masks = self._get_norge_i_bilder_polygon_masks(verbose=verbose)
+        for tile in self.tiles:
+            mask = masks.get(tile["name"], None)
+            tile["geometry"] = mask
+
+    def _get_norge_i_bilder_polygon_masks(self, verbose: bool):
+        from owslib.util import ServiceException
+        from owslib.wms import WebMapService
+        from PIL import Image
+
+        relevant_names: dict[str, str] = {x["name"]: x["bbox"] for x in self.tiles}
+        assert len(relevant_names), relevant_names
+
+        url = "https://wms.geonorge.no/skwms1/wms.nib-mosaikk?SERVICE=WMS&REQUEST=GetCapabilities"
+        wms = WebMapService(url, version="1.3.0")
+        out = {}
+        # ttiles = {wms[layer].title: [] for layer in list(wms.contents)}
+        # for layer in list(wms.contents):
+        #     if wms[layer].title not in relevant_names:
+        #         continue
+        #     ttiles[wms[layer].title].append(layer)
+        # import pandas as pd
+
+        # df = pd.Series(ttiles).to_frame("title")
+        # df["n"] = df["title"].str.len()
+        # df = df.sort_values("n")
+        # for x in df["title"]:
+        #     if len(x) == 1:
+        #         continue
+        #     bounds = {tuple(wms[layer].boundingBoxWGS84) for layer in x}
+        #     if len(bounds) <= 1:
+        #         continue
+        #     print()
+        #     for layer in x:
+        #         print(layer)
+        #         print(wms[layer].title)
+        #         bbox = wms[layer].boundingBoxWGS84
+        #         print(bbox)
+
+        for layer in list(wms.contents):
+            title = wms[layer].title
+            if title not in relevant_names:
+                continue
+            bbox = wms[layer].boundingBoxWGS84
+            bbox = tuple(to_gdf(bbox, crs=4326).to_crs(25832).total_bounds)
+
+            existing_bbox = relevant_names[title]
+            existing_bbox = to_gdf(existing_bbox, crs=4326).to_crs(25832).union_all()
+            if not to_shapely(bbox).intersects(existing_bbox):
+                continue
+            diffx = bbox[2] - bbox[0]
+            diffy = bbox[3] - bbox[1]
+            width = int(diffx / 40)
+            height = int(diffy / 40)
+            if not bbox:
+                continue
+            try:
+                img = wms.getmap(
+                    layers=[layer],
+                    styles=[""],  # Empty unless you know the style
+                    srs="EPSG:25832",
+                    bbox=bbox,
+                    size=(width, height),
+                    format="image/jpeg",
+                    transparent=True,
+                    bgcolor="#FFFFFF",
+                )
+            except (ServiceException, AttributeError) as e:
+                if verbose:
+                    print(type(e), e)
+                continue
+
+            arr = np.array(Image.open(BytesIO(img.read())))
+            if not np.sum(arr):
+                continue
+
+            band = Band(
+                np.where(np.any(arr != 0, axis=-1), 1, 0), bounds=bbox, crs=25832
+            )
+            polygon = band.to_geopandas()[lambda x: x["value"] == 1].geometry.values
+            polygon = make_valid(polygons(get_exterior_ring(polygon)))
+            polygon = make_valid(set_precision(polygon, 1))
+            polygon = make_valid(simplify(polygon, 100))
+            polygon = make_valid(set_precision(polygon, 1))
+            polygon = GeoSeries(polygon, crs=25832).to_crs(4326)
+            if verbose:
+                print(f"Layer name: {layer}")
+                print(f"Title: {wms[layer].title}")
+                print(f"Bounding box: {wms[layer].boundingBoxWGS84}")
+                print(f"polygon: {polygon}")
+                print("-" * 40)
+
+            for x in [0, 0.1, 0.001, 1]:
+                try:
+                    out[title] = make_valid(polygon.buffer(x).make_valid().union_all())
+                except GEOSException:
+                    pass
+                break
+
+        return out
+
+    def get_tiles(self, mask: Any, max_zoom: int = 40) -> list[folium.WmsTileLayer]:
+        """Get all Norge i bilder tiles intersecting with a mask (bbox or polygon)."""
         if self.tiles is None:
             self.load_tiles()
 
-        all_tiles = {}
-
-        bbox = to_shapely(bbox)
+        if not isinstance(mask, (GeoSeries | GeoDataFrame | Geometry)):
+            mask = to_shapely(mask)
 
         if isinstance(self.show, bool):
             show = self.show
         else:
             show = False
 
-        for tile in self.tiles:
-            if not tile["bbox"] or not tile["bbox"].intersects(bbox):
-                continue
-
-            name = tile["name"]
-
-            if (
-                not name
-                or not any(year in name for year in self.years)
-                or (
-                    self.contains
-                    and not any(re.search(x, name.lower()) for x in self.contains)
-                )
-                or (
-                    self.not_contains
-                    and any(re.search(x, name.lower()) for x in self.not_contains)
-                )
-            ):
-                continue
-
-            all_tiles[name] = folium.WmsTileLayer(
+        relevant_tiles = self._filter_tiles(mask)
+        tile_layers = {
+            name: folium.WmsTileLayer(
                 url="https://wms.geonorge.no/skwms1/wms.nib-prosjekter",
                 name=name,
                 layers=name,
@@ -161,16 +261,37 @@ def get_tiles(self, bbox: Any, max_zoom: int = 40) -> list[folium.WmsTileLayer]:
                 show=show,
                 max_zoom=max_zoom,
             )
+            for name in relevant_tiles["name"]
+        }
+
+        if not len(tile_layers):
+            return tile_layers
 
         if isinstance(self.show, int):
-            tile = all_tiles[list(all_tiles)[self.show]]
+            tile = tile_layers[list(tile_layers)[self.show]]
             tile.show = True
         elif isinstance(self.show, Iterable):
             for i in self.show:
-                tile = all_tiles[list(all_tiles)[i]]
+                tile = tile_layers[list(tile_layers)[i]]
                 tile.show = True
 
-        return all_tiles
+        return tile_layers
+
+    def _filter_tiles(self, mask):
+        """Filter relevant dates with pandas and geopandas because fast."""
+        df = pd.DataFrame(self.tiles)
+        filt = (df["name"].notna()) & (df["year"].str.contains("|".join(self.years)))
+        if self.contains:
+            for x in self.contains:
+                filt &= df["name"].str.contains(x)
+        if self.not_contains:
+            for x in self.not_contains:
+                filt &= ~df["name"].str.contains(x)
+        df = df[filt]
+        geoms = np.where(df["geometry"].notna(), df["geometry"], df["bbox"])
+        geoms = GeoSeries(geoms)
+        assert geoms.index.is_unique
+        return df.iloc[sfilter(geoms, mask).index]
 
     def __post_init__(self) -> None:
         """Fix typings."""
@@ -195,7 +316,11 @@ def __post_init__(self) -> None:
                 return
             self.tiles = [
                 {
-                    key: value if key != "bbox" else shapely.wkt.loads(value)
+                    key: (
+                        value
+                        if key not in ["bbox", "geometry"]
+                        else shapely.wkt.loads(value)
+                    )
                     for key, value in tile.items()
                 }
                 for tile in self.tiles

src/sgis/maps/wms.py

@@ -1,5 +1,4 @@
 # %%
-import datetime
 import json
 import os
 import warnings
@@ -117,10 +116,19 @@ def test_explore(points_oslo, roads_oslo):
     roads["geometry"] = roads.buffer(3)
 
     r300 = roads.clip(p.buffer(300))
-
     r200 = roads.clip(p.buffer(200))
     r100 = roads.clip(p.buffer(100))
 
+    e = sg.explore(
+        r300,
+        "meters",
+        r100,
+        bygdoy=7000,
+        wms=sg.NorgeIBilderWms(
+            years=sg.maps.wms.JSON_YEARS, not_contains="sentinel", show=-1
+        ),
+    )
+    assert isinstance(next(iter(e.wms)), sg.NorgeIBilderWms)
     e = sg.explore(
         r300,
         "meters",
@@ -285,10 +293,13 @@ def not_test_wms_json():
     print(
         "IMPORTANT: if you run this function, make sure to change the global variable JSON_YEARS in wms.py"
     )
-    wms = sg.NorgeIBilderWms(
-        years=range(1999, datetime.datetime.now().year + 1), _use_json=False
-    )
-    wms.load_tiles()
+    years = sg.maps.wms.JSON_YEARS
+    wms = sg.NorgeIBilderWms(years=years, _use_json=False)
+    wms.load_tiles(verbose=True)
+
+    with open(sg.maps.wms.JSON_PATH) as file:
+        print(pd.DataFrame(file).year.value_counts())
+
     try:
         os.remove(sg.maps.wms.JSON_PATH)
     except FileNotFoundError:
@@ -297,7 +308,7 @@ def not_test_wms_json():
         json.dump(
             [
                 {
-                    key: value if key != "bbox" else value.wkt
+                    key: value if key not in ["bbox", "geometry"] else value.wkt
                     for key, value in tile.items()
                 }
                 for tile in wms.tiles
@@ -306,16 +317,23 @@ def not_test_wms_json():
             ensure_ascii=False,
         )
 
+    wms = sg.NorgeIBilderWms(years=years)
+    for tile in wms.tiles:
+        bbox = sg.to_shapely(tile["bbox"])
+        polygon = sg.to_shapely(tile["geometry"])
+        assert bbox.intersects(polygon)
+
 
 def main():
 
     from oslo import points_oslo
     from oslo import roads_oslo
 
-    not_test_wms_json()
+    # not_test_wms_json()
 
     test_explore(points_oslo(), roads_oslo())
     test_image_collection()
+
     # not_test_explore(points_oslo(), roads_oslo())