wrpylib/wrgeojson.py

   1 """Winterrodeln maps can be represented in the GeoJSON format - this is referred to as wrgeojson.
   2 This module contains related helper functions.
   3
   4 Recommended/used Python packages
   5
   6 * To work with a wrgeojson map, it is recommended but not necessary
   7   to use the Python [geojson](https://github.com/jazzband/geojson) package.
   8   This is especially useful as it provides type annotations for the GeoJSON types.
   9 * For coordinate transformation, the [pyproj](https://pyproj4.github.io/pyproj/stable/) package is recommended.
  10 * For math with geometric shapes, [shapely](https://github.com/shapely/shapely) is recommended.
  11
  12
  13 The definition is kept in sync with the definition of wrmap, available at
  14 git+ssh://gitolite@git.toastfreeware.priv.at/philipp/winterrodeln/mediawiki_extensions/wrmap.git
  15 (in wrmap.body.php)
  16
  17 A Winterrodeln  map is represented as `feature collection`.
  18 Feature collections must not have a JSON name "properties",
  19 but accidentally this name was used to store map related Winterrodeln specific
  20 properties.
  21 Those properties were:
  22
  23 * 'lat': float: latitude of map-center, optional.
  24 * 'lon': float: longitude of map-center, optional.
  25 * 'zoom': int: zoom level of the map (google zoom levels). optional.
  26 * 'width': int: width of the map in pixel. optional (100% if omitted)
  27 * 'height': int: height of the map in pixel. optional.
  28 * 'under_construction': bool: whether the map is unfinished, optional.
  29 * 'simplified': bool: whether the map was simplified (usually with simplify_ways()), optional.
  30
  31 The new maps will stick to the GeoJSON definition https://datatracker.ietf.org/doc/html/rfc7946
  32 and use an optional foreign member object with the name "wr_properties", where the following optional names
  33 are allowed:
  34
  35 * 'lat': float: latitude of map-center, optional and deprecated in favor of vbox.
  36 * 'lon': float: longitude of map-center, optional and deprecated in favor of vbox.
  37 * 'zoom': int: zoom level of the map (google zoom levels), optional and deprecated in favor of vbox.
  38 * 'width': int: width of the map in pixel, optional and deprecated in favor of vbox (100% if omitted).
  39 * 'height': int: height of the map in pixel, optional and deprecated in favor of vbox.
  40 * 'vbox': array of 4 floats, optional. Same meaning as bbox but contains the map region that should be shown by default.
  41 * 'under_construction': bool: whether the map is unfinished, optional. Default: False
  42 * 'simplified': bool: whether the map was simplified (usually with simplify_ways()), optional.
  43   Currently all maps are converted to be not simplified anymore so in the future this name will be deprecated.
  44
  45
  46 For each feature inside the feature collection the following rules apply:
  47
  48 * The supported feature geometries are of type Point or LineString.
  49 * Each feature has a property called type with a string value.
  50 * The type of each Point feature is one of
  51   * 'gasthaus'
  52   * 'haltestelle'
  53   * 'parkplatz'
  54   * 'achtung'
  55   * 'foto'
  56   * 'verleih'
  57   * 'punkt'
  58
  59 * The type of each LineString feature is one of
  60   * 'rodelbahn'
  61   * 'alternative'
  62   * 'gehweg'
  63   * 'lift'
  64   * 'anfahrt'
  65   * 'linie'
  66
  67
  68 Example:
  69
  70 {
  71     "type": "FeatureCollection",
  72     "features": [
  73         {
  74             "type": "Feature",
  75             "geometry": {"type": "Point", "coordinates": [11.42134, 47.302911]},
  76             "properties": {
  77                 "type": "gasthaus",
  78                 "name": "Rumer Alm",
  79                 "wiki": "Rumer Alm (Gasthaus)"
  80             }
  81         },
  82         {
  83             "type": "Feature",
  84             "geometry": {"type": "Point", "coordinates": [11.448077, 47.288761]},
  85             "properties": {
  86                 "type": "haltestelle",
  87                 "name": "Rum Sanatorium"
  88             }
  89         },
  90         {
  91             "type": "Feature",
  92             "geometry": {
  93                 "type": "LineString",
  94                 "coordinates": [
  95                     [11.435297, 47.297118],
  96                     [11.434738, 47.297678],
  97                     [11.434973, 47.299145]
  98                 ]
  99             },
 100             "properties": {
 101                 "type": "gehweg"
 102             }
 103         },
 104         {
 105             "type": "Feature",
 106             "geometry": {
 107                 "type": "LineString",
 108                 "coordinates": [
 109                     [11.420065, 47.29844],
 110                     [11.418975, 47.298389],
 111                     [11.417893, 47.29784]
 112                 ]
 113             },
 114             "properties": {
 115                 "type": "alternative"
 116             }
 117         },
 118         {
 119             "type": "Feature",
 120             "geometry": {
 121                 "type": "LineString",
 122                 "coordinates": [
 123                     [11.421153, 47.302928],
 124                     [11.418656, 47.303768],
 125                     [11.439608, 47.293187],
 126                     [11.441071, 47.291241]
 127                 ]
 128             },
 129             "properties": {
 130                 "type": "rodelbahn"
 131             }
 132         }
 133     ],
 134     "wr_properties": {
 135         "lat": 47.29626785,
 136         "lon": 11.4327404,
 137         "zoom": 14,
 138         "simplified": True
 139     }
 140 }
 141
 142
 143 Recipies for working with GeoJSON:
 144
 145 Load GeoJSON:
 146
 147 ```
 148 import geojson
 149 with open('wrgeojson.geojson') as fp:
 150     gj = geojson.
 151 ```
 152
 153 Distance between two points:
 154 ```
 155 from pyproj import Geod # https://pyproj4.github.io/pyproj/stable/
 156 from shapely.geometry import Point, LineString
 157
 158 geod = Geod(ellps="WGS84")
 159 p1 = Point(11.404997, 47.283445)
 160 p2 = Point(11.405752, 47.284616)
 161 line_string = LineString([p1, p2])
 162 ```
 163
 164
 165 Local CRS (meters)
 166 ```
 167 from pyproj import CRS, Transformer
 168 crs_from = CRS.from_proj4("+proj=latlon")
 169 crs_to = CRS.from_proj4('+proj=merc +lat_ts=47.3')
 170 transformer = Transformer.from_crs(crs_from, crs_to, always_xy=True)
 171 p1t = transformer.transform(*list(*p1.coords))
 172 p2t = transformer.transform(*list(*p2.coords))
 173 import math
 174 math.sqrt((p2t[0]-p1t[0])**2 + (p2t[1]-p1t[1])**2)  # 142.12323235553566
 175 ```
 176
 177 Work with shapely:
 178 ```
 179 p1s = shapely.geometry.shape(p1)
 180 print(p1s.wkt)
 181 ```
 182
 183 Distance/similarity between two lines:
 184 shapely.hausdorff_distance might be used.
 185 """
 186 import math
 187 from collections import Counter
 188 from copy import deepcopy
 189 from enum import Enum
 190 from typing import Tuple, List
 191
 192 import geojson.utils
 193 import shapely
 194 import shapely.geometry
 195 from geojson import FeatureCollection, Feature, Point, GeoJSON
 196 from pyproj import CRS, Transformer
 197
 198
 199 class WrPointFeatureType(Enum):
 200     gastronomy = 'gasthaus'
 201     bus_stop = 'haltestelle'
 202     parking = 'parkplatz'
 203     warning = 'achtung'
 204     foto = 'foto'
 205     rental = 'verleih'
 206     point = 'punkt'
 207
 208
 209 class WrLineStringFeatureType(Enum):
 210     sled = 'rodelbahn'
 211     alternative = 'alternative'
 212     walk = 'gehweg'
 213     lift = 'lift'
 214     drive = 'anfahrt'
 215     line = 'linie'
 216
 217
 218 def line_point_dist(lon_lat_1: Tuple[float, ...], lon_lat_2: Tuple[float, ...],
 219                     lon_lat: Tuple[float, ...]) -> float:
 220     """Returns the distance in meter between the line defined by
 221     lon_lat_1, lon_lat_2
 222     and the point defined by lon_lat."""
 223     # Do the calculation in m
 224     earth_radius = 6367500.  # earth radius in m
 225     deg_to_rad = math.pi / 180.
 226     lat_to_meter = deg_to_rad * earth_radius
 227
 228     if lon_lat_1 == lon_lat_2:
 229         dx = (lon_lat[0] - lon_lat_1[0]) * math.cos(lon_lat_1[0]) * lat_to_meter
 230         dy = (lon_lat[1] - lon_lat_1[1]) * lat_to_meter
 231         return math.sqrt(dx ** 2 + dy ** 2)
 232
 233     # d = lon_lat_2 - lon_lat_1
 234     # e = lon_lat - lon_lat_1
 235
 236     mean_lat_rad = (lon_lat_1[1] + lon_lat_2[1] + lon_lat[1]) / 3. * deg_to_rad
 237
 238     dx = (lon_lat_2[0] - lon_lat_1[0]) * math.cos(mean_lat_rad) * lat_to_meter
 239     dy = (lon_lat_2[1] - lon_lat_1[1]) * lat_to_meter
 240     ex = (lon_lat[0] - lon_lat_1[0]) * math.cos(mean_lat_rad) * lat_to_meter
 241     ey = (lon_lat[1] - lon_lat_1[1]) * lat_to_meter
 242
 243     ds = dx ** 2 + dy ** 2
 244     es = ex ** 2 + ey ** 2
 245     xs = es - (dx * ex + dy * ey) ** 2 / ds
 246     if xs < 0:
 247         return 0.  # due to rounding errors xs might be e.g. -1e-23
 248     return math.sqrt(xs)
 249
 250
 251 def simplify_way(lon_lat_list: List[Tuple[float, ...]], max_dist: float):
 252     """Simplify lon_lat_list by removing points in-place. The introduced error must not exceed max_dist (meter)."""
 253     dist_list = []  # dist_list[i] is the distance between
 254     # the line lon_lat_list[i] <-> lon_lat_list[i+2]
 255     # and the point lon_lat_list[i+1].
 256     # len(dist_list) == len(lon_lat_list)-2
 257
 258     # Initialize dist_list
 259     for i in range(len(lon_lat_list) - 2):
 260         dist_list.append(line_point_dist(lon_lat_list[i], lon_lat_list[i + 2], lon_lat_list[i + 1]))
 261
 262     while len(dist_list) > 0:
 263         min_index = min(range(len(dist_list)), key=dist_list.__getitem__)  # same as index(min(dist_list))
 264         if dist_list[min_index] > max_dist:
 265             return
 266         del lon_lat_list[min_index + 1]
 267         del dist_list[min_index]
 268
 269         if min_index > 0:
 270             dist_list[min_index - 1] = line_point_dist(lon_lat_list[min_index - 1], lon_lat_list[min_index + 1],
 271                                                        lon_lat_list[min_index],)
 272         if min_index < len(dist_list):
 273             dist_list[min_index] = line_point_dist(lon_lat_list[min_index], lon_lat_list[min_index + 2],
 274                                                    lon_lat_list[min_index + 1])
 275
 276
 277 DEFAULT_MAX_DIST_ERROR_M = 15.
 278
 279
 280 def simplify_ways(wrgeojson: FeatureCollection, max_dist: float = DEFAULT_MAX_DIST_ERROR_M):
 281     ways = [feature for feature in wrgeojson['features'] if feature["geometry"]["type"] == "LineString"]
 282     for way in ways:
 283         simplify_way(way['geometry']['coordinates'], max_dist)
 284
 285
 286 def get_geometry_center(wrgeojson: GeoJSON) -> Point:
 287     return Point(get_geometry_envelope(wrgeojson).centroid.coords[0])
 288
 289
 290 def get_geometry_points(wrgeojson: GeoJSON) -> shapely.geometry.LineString:
 291     return shapely.geometry.LineString(c[:2] for c in geojson.coords(wrgeojson))
 292
 293
 294 def get_geometry_envelope(wrgeojson: GeoJSON) -> shapely.geometry.Polygon:
 295     return get_geometry_points(wrgeojson).envelope
 296
 297
 298 def get_geometry_bbox(wrgeojson: GeoJSON) -> Tuple[float, float, float, float]:
 299     return get_geometry_envelope(wrgeojson).bounds
 300
 301
 302 def sort_features(wrgeojson: FeatureCollection) -> FeatureCollection:
 303     """Sort the features of the sledruns in preferred order of rendering"""
 304     def sort_key(feature: Feature) -> Tuple[int, int]:
 305         geometry_type = feature['geometry']['type']
 306         feature_type = feature['properties'].get('type')
 307         if geometry_type == 'LineString':
 308             # left types are drawn first
 309             order = [None, 'linie', 'anfahrt', 'gehweg', 'alternative', 'rodelbahn', 'lift']
 310             return 1, order.index(feature_type)
 311         else:
 312             order = ['gasthaus', 'haltestelle', 'parkplatz', 'warnung', 'punkt']  # the left types are coming first
 313             return 0, order.index(feature_type)
 314     wrgeojson.features = sorted(wrgeojson.features, key=sort_key)
 315     return wrgeojson
 316
 317
 318 def join_wrgeojson_ways(wrgeojson: FeatureCollection, preserve_crossing=False):
 319     """Unite ways of the same type of sled runs."""
 320     # Get crossing (lon, lat) tuples as list.
 321     # Crossings are positions (lon, lat) that occur at least three way-edges.
 322     ways = [feature for feature in wrgeojson['features'] if feature["geometry"]["type"] == "LineString"]
 323     coord_count = Counter(geojson.utils.coords(ways))
 324     crossing_list = [coord for coord, count in coord_count.items() if count >= 3]
 325
 326     i = 0  # index within ways
 327     while i < len(ways):
 328         # search for ways to combine (they have to have the same type)
 329         if len(ways[i]['geometry']['coordinates']) == 0:
 330             i += 1
 331             continue
 332         j = i + 1
 333         while j < len(ways):
 334             # comparing i and j
 335             if ways[i]['properties']['type'] != ways[j]['properties']['type'] \
 336                     or len(ways[j]['geometry']['coordinates']) == 0:
 337                 j += 1
 338                 continue
 339             lon_lat = ways[i]['geometry']['coordinates'][-1]
 340             if not preserve_crossing or lon_lat not in crossing_list:
 341                 del_j = False
 342                 if lon_lat == ways[j]['geometry']['coordinates'][0]:
 343                     # Append way
 344                     ways[i]['geometry']['coordinates'].extend(ways[j]['geometry']['coordinates'][1:])
 345                     del_j = True
 346                 elif lon_lat == ways[j]['geometry']['coordinates'][-1]:
 347                     # Append reversed way
 348                     ways[i]['geometry']['coordinates'].extend(ways[j]['geometry']['coordinates'][-2::-1])
 349                     del_j = True
 350                 else:
 351                     lon_lat = ways[i]['geometry']['coordinates'][0]
 352                     if lon_lat == ways[j]['geometry']['coordinates'][-1]:
 353                         # Prepend way
 354                         ways[i]['geometry']['coordinates'] = ways[j]['geometry']['coordinates'][:-1] + \
 355                                                              ways[i]['geometry']['coordinates']
 356                         del_j = True
 357                     elif lon_lat == ways[j]['geometry']['coordinates'][0]:
 358                         # Prepend reversed way
 359                         ways[i]['geometry']['coordinates'] = ways[j]['geometry']['coordinates'][:0:-1] + \
 360                                                              ways[i]['geometry']['coordinates']
 361                         del_j = True
 362                 if del_j:
 363                     ways[j]['geometry']['coordinates'].clear()
 364                 j += 1
 365         i += 1
 366     wrgeojson['features'] = [f for f in wrgeojson['features']
 367                              if f['geometry']['type'] != 'LineString' or len(f['geometry']['coordinates']) > 0]
 368
 369
 370 class LanToToMetricTransformer:
 371     def __init__(self, latitude: float):
 372         crs_from = CRS.from_proj4("+proj=latlon")
 373         crs_to = CRS.from_proj4(f'+proj=merc +lat_ts={latitude:.3f}')
 374         self.transformer = Transformer.from_crs(crs_from, crs_to, always_xy=True)
 375
 376     def tuple_transformer(self, lon_lat: Tuple[float, float]) -> Tuple[float, float]:
 377         return self.transformer.transform(lon_lat[0], lon_lat[1])
 378
 379     def geojson_lon_lat_to_metric(self, gj: GeoJSON) ->GeoJSON:
 380         gj = deepcopy(gj)
 381         geojson.utils.map_tuples(self.tuple_transformer, gj)
 382         return gj
 383
 384
 385 def geojson_lon_lat_to_metric(gj: GeoJSON) -> GeoJSON:
 386     """Converts the coordinates of GeoJSON FeatureCollection as well as other GeoJSON objects to a metric
 387     coordinate system."""
 388     center = get_geometry_center(gj)
 389     transformer = LanToToMetricTransformer(center.coordinates[1])
 390     return transformer.geojson_lon_lat_to_metric(gj)