Tobler

osm-elevation.py

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import geopyspark as gps\n",
    "import numpy as np\n",
    "from pyspark import SparkContext, StorageLevel\n",
    "from pyspark.sql import SparkSession\n",
    "from geonotebook.wrappers import TMSRasterData\n",
    "\n",
    "conf = gps.geopyspark_conf(appName=\"gps-osm-ingest\", master='local[8]')\n",
    "conf.set(\"spark.hadoop.yarn.timeline-service.enabled\", False)\n",
    "conf.set('spark.ui.enabled', True)\n",
    "conf.set('spark.default.parallelism', 8)\n",
    "conf.set('spark.master.memory', '9500M')\n",
    "\n",
    "sc = SparkContext(conf=conf)\n",
    "pysc = gps.get_spark_context()\n",
    "session = SparkSession.builder.config(conf=pysc.getConf()).enableHiveSupport().getOrCreate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "file_uri = \"file:/Users/eugene/Downloads/isle-of-man.orc\"\n",
    "osm_dataframe = session.read.orc(file_uri)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "osm = gps.vector_pipe.osm_reader.from_dataframe(osm_dataframe)\n",
    "lines = osm.get_line_features_rdd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Feature(geometry=<shapely.geometry.linestring.LineString object at 0x10e768438>, properties=Properties(element_id=112582296, user='FireCred', uid=594323, changeset=11161182, version=1, minor_version=1, timestamp=datetime.datetime(2012, 3, 30, 21, 15, 59, tzinfo=tzutc()), visible=True, tags={'service': 'alley', 'highway': 'service'}))]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lines.take(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "highways = lines.filter(\n",
    "    lambda feature: 'highway' in feature.properties.tags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Feature(geometry=<shapely.geometry.linestring.LineString object at 0x10e768940>, properties=Properties(element_id=112582296, user='FireCred', uid=594323, changeset=11161182, version=1, minor_version=1, timestamp=datetime.datetime(2012, 3, 30, 21, 15, 59, tzinfo=tzutc()), visible=True, tags={'service': 'alley', 'highway': 'service'})),\n",
       " Feature(geometry=<shapely.geometry.linestring.LineString object at 0x10e768908>, properties=Properties(element_id=112582296, user='manksman', uid=364604, changeset=8102929, version=1, minor_version=0, timestamp=datetime.datetime(2011, 5, 10, 12, 29, 2, tzinfo=tzutc()), visible=True, tags={'service': 'alley', 'highway': 'service'}))]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "highways.take(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['50 mph',\n",
       " '10 mph',\n",
       " '5 mph',\n",
       " '20',\n",
       " '50',\n",
       " '30 mph',\n",
       " '60 mph',\n",
       " '3mph',\n",
       " '40 mph',\n",
       " '15 mph',\n",
       " 'unrestricted',\n",
       " '5mpg',\n",
       " '5',\n",
       " 'national',\n",
       " '20 mph',\n",
       " '15mph',\n",
       " 'none',\n",
       " '50mph',\n",
       " '30mph',\n",
       " '60mph',\n",
       " '32',\n",
       " '20mph',\n",
       " '40mph',\n",
       " '10mph',\n",
       " '5mph',\n",
       " '3 mph']"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "highways.filter(lambda f: 'maxspeed' in f.properties.tags).\\\n",
    "    map(lambda f: f.properties.tags['maxspeed']).distinct().take(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "highway_features = highways.map(\n",
    "    lambda feature: gps.Feature(feature.geometry, gps.CellValue(2, 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "highway_raster = gps.geotrellis.rasterize_features(\n",
    "    features = highway_features,\n",
    "    crs = 4326,\n",
    "    zoom = 10,\n",
    "    cell_type=gps.CellType.INT8RAW,\n",
    "    num_partitions = 16\n",
    ").persist(StorageLevel.MEMORY_AND_DISK)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# insert block that renders the lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<promise.promise.Promise at 0x10de67ba8>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "color_map = gps.ColorMap.from_colors(\n",
    "    breaks = np.arange(0, 3, 1), \n",
    "    color_list = gps.get_colors_from_matplotlib('magma'))\n",
    "\n",
    "osm_wm = highway_raster.tile_to_layout(\n",
    "    gps.GlobalLayout(), target_crs=3857)\n",
    "layer = gps.TMS.build(osm_wm.pyramid(), color_map)\n",
    "\n",
    "M.add_layer(TMSRasterData(layer), name=\"OSM\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "isle_extent = gps.Extent(xmin=-4.921879863281248, ymin=53.964837753906245, \n",
    "                         xmax=-4.218754882812476, ymax=54.4921814453125)\n",
    "#isle_extent = highway_raster.layer_metadata.extent\n",
    "\n",
    "# Read SRTM for the same extent as rasterized OSM features\n",
    "srtm = gps.query(\n",
    "    uri=\"s3://geotrellis-test/dg-srtm\",\n",
    "    layer_name=\"srtm-wsg84-gps\", \n",
    "    layer_zoom = 0,\n",
    "    query_geom = isle_extent,\n",
    "    query_proj = 4326,\n",
    "    num_partitions = 16\n",
    ")\n",
    "\n",
    "# Tile SRTM layer to same layout as rasterized OSM features\n",
    "tiled_srtm = srtm.tile_to_layout(highway_raster.layer_metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compensate for cell size at 55N\n",
    "z_factor = 0.000015935\n",
    "\n",
    "tobler_raster = tiled_srtm.\\\n",
    "    focal(gps.Operation.SLOPE, param_1=z_factor).\\\n",
    "    tobler()\n",
    "    \n",
    "friction = tobler_raster - highway_raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reprojected = friction.tile_to_layout(\n",
    "    target_crs = 3857,\n",
    "    layout = gps.GlobalLayout(tile_size=256),\n",
    "    resample_method = gps.ResampleMethod.MAX\n",
    ").persist(StorageLevel.MEMORY_AND_DISK_SER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pyramid = reprojected.pyramid().persist(StorageLevel.MEMORY_AND_DISK_SER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<promise.promise.Promise at 0x10dea68d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "color_map = gps.ColorMap.from_colors(\n",
    "    breaks = np.arange(0, 7.5, 0.2), \n",
    "    color_list = gps.get_colors_from_matplotlib('magma'))\n",
    "\n",
    "layer = gps.TMS.build(pyramid, color_map)\n",
    "\n",
    "M.add_layer(TMSRasterData(layer), name=\"ToblerOSM\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.8699236113597981, 5.036742124615245)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyramid.get_histogram().min_max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<promise.promise.Promise at 0x10de90b38>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.remove_layer(M.layers[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TiledRasterLayer(layer_type=LayerType.SPATIAL, zoom_level=13, is_floating_point_layer=False)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reprojected.unpersist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Write pyramid to the catalog\n",
    "for layer in pyrmaid.levels.values():\n",
    "    gps.write(\"s3://geotrellis-test/dg-osm-test\", \"gps-osm-ingest\", layer)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Geonotebook (Python 3)",
   "language": "python",
   "name": "geonotebook3"
  },
  "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}