ShapesFilter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare DAV Wildschutzgebiete (epsg:4326 shape file) with mapped Wildschutzgebiete (XCTrails geojson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\r\n",
    "from shapely.geometry import shape, mapping\r\n",
    "from shapely.ops import transform\r\n",
    "from shapely.geometry import Polygon\r\n",
    "\r\n",
    "import copy\r\n",
    "\r\n",
    "cat = 'Wildschutzgebiet'\r\n",
    "\r\n",
    "mappedFile = 'E:/OSM/Schutzgebiete/Schongebiete.geojson' # epsg:4326\r\n",
    "\r\n",
    "shapeFile = 'E:/OSM/Schutzgebiete/200126_Schutzgebiete_By-Karten/Schutzgebiete_BY-Karten-epsg4326.shp' # epsg:4326\r\n",
    "outFile = f'E:/OSM/Schutzgebiete/200126_Schutzgebiete_By-Karten/Schutzgebiete_BY-Karten-{cat}.shp' # epsg:4326"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OSM-Gebiete in E:/OSM/Schutzgebiete/Schongebiete.geojson: 557\n"
     ]
    }
   ],
   "source": [
    "mapped_features = []\r\n",
    "with fiona.open(mappedFile) as input:\r\n",
    "    oldCrs = input.crs\r\n",
    "    for feat in input:\r\n",
    "        if feat['geometry'] != None and len(feat['geometry']['coordinates'][0]) > 2:\r\n",
    "            mapped_features.append(feat)\r\n",
    "size = len(mapped_features)\r\n",
    "print(f\"OSM-Gebiete in {mappedFile}: {size}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find a given poly in a list of geometry features and return its (mapped) properties\r\n",
    "def find_poly(poly, features):\r\n",
    "    poly = poly.buffer(0.0001)\r\n",
    "    for f in features:\r\n",
    "        p = Polygon(shape(f['geometry'])).buffer(0.0001)\r\n",
    "        iou = poly.intersection(p).area / poly.union(p).area\r\n",
    "        if iou > 0.995:\r\n",
    "            return f['properties']\r\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\"Wirtsalm\" Betretungsverbot: 01.11. bis 31.03.:\n",
      "\tway/923656164 type2 Wildschutzgebiet Wirtsalm Betretungsverbot: 01.11. bis 31.03. no @ (Nov - Mar)\n",
      "\"Feichteck\" Betretungsverbot: 15.10. bis 15.05.:\n",
      "\tway/925219223 type2 Wildschutzgebiet Feichteck Betretungsverbot: 15.10. bis 15.05. no @ (Oct 15 - May 15)\n",
      "\"Grattenbach\" Betretungsverbot: 01.11. bis 31.05.:\n",
      "\tway/925194626 type2 Wildschutzgebiet Grattenbach Betretungsverbot: 01.11. bis 31.05. no @ (Nov - May)\n",
      "\"Gründberg\" Betretungsverbot: 01.11. bis 30.04.:\n",
      "\tway/925185305 type2 Wildschutzgebiet Gründberg Betretungsverbot: 01.11. bis 30.04. no @ (Nov 1 - Apr 30)\n",
      "\"Bärenfalle\" Betretungsverbot: 15.01. bis 30.04.:\n",
      "\tway/923832592  Wald-Wild-Schongebiet Bärenfalle Wald-Wild-Schongebiet, mit der Bitte dieses nicht zu betreten oder zu befahren. None\n",
      "\"Halmrauth\" Betretungsverbot: 01.11. bis 15.04.:\n",
      "\tway/925199679 type2 Wildschutzgebiet Halmrauth Betretungsverbot: 01.11. bis 15.04. no @ (Nov 1 - Apr 15)\n",
      "\"Bärenleiten\" Betretungsverbot: 01.11. bis 15.04.:\n",
      "\tway/925199677  Wald-Wild-Schongebiet Skitouren- und Schneeschuhgeher: Bitte nicht betreten oder befahren! None\n",
      "\"Schneckenwurf\" Betretungsverbot: 01.11. bis 15.04.:\n",
      "\tway/925219220 type2 Wildschutzgebiet Schneckenwurf Betretungsverbot: 01.11. bis 15.04. no @ (Nov 01-Apr 15)\n",
      "\"Laingraben\" Betretungsverbot: 01.11. bis 15.04.:\n",
      "\tway/925219222 type2 Wildschutzgebiet Laingraben Betretungsverbot: 01.11. bis 15.04. no @ (Nov 1-Apr 15)\n",
      "\"Hachel\" Betretungsverbot: 01.11. bis 15.04.:\n",
      "\tway/925198058 type2 Wildschutzgebiet Hachel Betretungsverbot: 01.11. bis 15.04. no @ (Nov 1 - Apr 15)\n",
      "\"Rotenkopf Süd\" Betretungsverbot: 15.11. bis 15.04.:\n",
      "\tway/925219218 type2 Wildschutzgebiet Rotenkopf Süd Betretungsverbot: 15.11. bis 15.04. no @ (Nov 15-Apr 15)\n",
      "\"Reintal\" Betretungsverbot: 15.11. bis 15.04.:\n",
      "\tway/925219219 type2 Wildschutzgebiet Reintal Betretungsverbot: 15.11. bis 15.04. no @ (Nov 15-Apr 15)\n",
      "\"Siebenhütten\" Betretungsverbot  01.12. bis 31.03., 12 Uhr bis 08 Uhr:\n",
      "\tway/923695934 type2 Wildschutzgebiet Hohlenstein Wildschutzgebiet mit Betretungsverbot 1.12. – 31.3. von 12:00 Uhr bis 8:00 Uhr (Gebiet darf nur von 8-12 Uhr betreten werden) no @ (Dec-Mar 12:00-08:00)\n",
      "\"Gescheuerwand\" Betretungsverbot: 01.12. bis 31.03.:\n",
      "\tway/925189143 type2 WildSchutzgebier Gescheuerwand Betretungsverbot: 01.12. bis 31.03. no @ (Dec - Mar)\n",
      "\"Wimbach\" Betretungsverbot: 15.10. bis 15.05.:\n",
      "\tway/925190292 type2 Wildschutzgebiet Wimbach Betretungsverbot: 15.10. bis 15.05. no @ (Oct 15 - May 15)\n",
      "\"Schwarzenkopf\" Betretungsverbot: ganzjährig:\n",
      "\tway/923695940 type1 Wald-Wild-Schongebiet Wildschutzgebiet mit ganzjährigem Betretungsverbot None\n",
      "\"Bleckstein Winterstube\" Betretungsverbot: 15.11. bis 31.05.:\n",
      "\tway/923695941 type2 Bleckstein Winterstube Wildschutzgebiet mit Betretungsverbot 15.11. – 31.5. no @ (Nov 15 - May 31)\n",
      "Gebiete in E:/OSM/Schutzgebiete/200126_Schutzgebiete_By-Karten/Schutzgebiete_BY-Karten-epsg4326.shp: 39\n"
     ]
    }
   ],
   "source": [
    "# read and 3D to 2D convert DAV shapefile\r\n",
    "newFeatures = []\r\n",
    "with fiona.open(shapeFile) as input:\r\n",
    "    schema = input.schema\r\n",
    "    crs = input.crs\r\n",
    "    driver = input.driver\r\n",
    "    for feat in input:\r\n",
    "        if 'Kategorie' in feat['properties'] and feat['properties']['Kategorie'] != None and feat['properties']['Kategorie'] == cat:\r\n",
    "            if feat['geometry'] != None:\r\n",
    "                if len(feat['geometry']['coordinates']) > 1:\r\n",
    "                    # multipolygons - code to be improved...\r\n",
    "                    for pfeat in feat['geometry']['coordinates']:\r\n",
    "                        try: # some are len(1) lists, some aren't ?!\r\n",
    "                            poly = Polygon(pfeat[0])\r\n",
    "                        except:\r\n",
    "                            poly = Polygon(pfeat)\r\n",
    "                        poly = transform(lambda x, y, z=None: (x, y), poly)\r\n",
    "                        feat2 = copy.deepcopy(feat)\r\n",
    "                        feat2['geometry'] = mapping(poly)\r\n",
    "                        newFeatures.append(feat2)\r\n",
    "                    continue\r\n",
    "                if len(feat['geometry']['coordinates'][0]) < 3:\r\n",
    "                    print(\"Skipping 2-point line\")\r\n",
    "                    continue\r\n",
    "                # transform 3D to 2D\r\n",
    "                poly = shape(feat['geometry'])\r\n",
    "                poly = transform(lambda x, y, z=None: (x, y), poly)\r\n",
    "                feat['geometry'] = mapping(poly)\r\n",
    "                mapped_props = find_poly(poly, mapped_features)\r\n",
    "                if not mapped_props:\r\n",
    "                    newFeatures.append(feat)\r\n",
    "                else:\r\n",
    "                    # print(mapped_props)\r\n",
    "                    # visual check in output:\r\n",
    "                    # - type1 oder type2\r\n",
    "                    # - name = Wildschutzgebiet ...\r\n",
    "                    # - Regelung passt mit description und access:conditional zusammen\r\n",
    "                    # WhiteList:\r\n",
    "                    # - Wald-Wild-Schongebiet Bärenfalle \r\n",
    "                    print(f\"{feat['properties']['Name']} {feat['properties']['Regelung']}:\")\r\n",
    "                    print(f\"\\t{mapped_props['@id']} {mapped_props['classification']} {mapped_props['name']} {mapped_props['description']} {mapped_props['access:conditional']}\")\r\n",
    "\r\n",
    "with fiona.open(outFile, 'w', crs={'init':'epsg:4326'}, driver='ESRI Shapefile', schema=schema) as out:\r\n",
    "    for f in newFeatures:\r\n",
    "        out.write(f)\r\n",
    "\r\n",
    "size = len(newFeatures);          \r\n",
    "print(f\"Gebiete in {shapeFile}: {size}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.2 64-bit ('.env': venv)",
   "metadata": {
    "interpreter": {
     "hash": "c1583a91ba8d9042e6938cef4a38a2659d1132ba00553427a7e20d2e070e8f70"
    }
   },
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}