prep_static_features.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import arcpy\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create GDB\n",
    "All data is remote, but some will be processed and saved to a file GDB. Create this GDB is it doesn't exist and set this as the workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if not os.path.exists('utah.gdb'):\n",
    "    arcpy.management.CreateFileGDB('.','utah.gdb')\n",
    "\n",
    "# workspace\n",
    "arcpy.env.workspace = r'./utah.gdb'\n",
    "arcpy.env.overwriteOutput = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define data URLs/Paths\n",
    "The below cell contains all the data needed. All data is remote except **collisions_csv** and **weather_stations_csv**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# AADT, Intersections, Billboards\n",
    "aadt_url = 'https://services1.arcgis.com/MdyCMZnX1raZ7TS3/arcgis/rest/services/Utah_AADT_2013/FeatureServer/0'\n",
    "intersections_url = 'https://maps.udot.utah.gov/arcgis/rest/services/FI_Mandli2012/MapServer/5'\n",
    "billboards_url = 'https://maps.udot.utah.gov/arcgis/rest/services/FI_Mandli2012/MapServer/2'\n",
    "\n",
    "collisions_csv = r'Utah crashes_2010-2017.csv'\n",
    "weather_stations_csv = 'utah_stations.csv'\n",
    "\n",
    "# Roads data\n",
    "major_roads = 'http://services1.arcgis.com/99lidPhWCzftIe9K/arcgis/rest/services/Utah_Major_Road_Centerlines_(Statewide)/FeatureServer/0'\n",
    "minor_roads = 'http://services1.arcgis.com/99lidPhWCzftIe9K/arcgis/rest/services/Utah_SGID_Road_Centerlines_(Statewide)/FeatureServer/0'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Copy the remote major roads features into local GDB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Copy the major road centerlines to the gdb\n",
    "_ = arcpy.management.CopyFeatures(in_features=major_roads,out_feature_class='centerlines')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Join AADT\n",
    "_ = arcpy.analysis.SpatialJoin('centerlines',\n",
    "                           aadt_url,\n",
    "                           'centerlines_aadt',\n",
    "                           match_option='closest',\n",
    "                           search_radius='200 feet')\n",
    "\n",
    "                           \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_ = arcpy.management.DeleteField('centerlines_aadt',['route_id','Shape_Leng','OBJECTID_12'])\n",
    "_ = arcpy.management.AlterField('centerlines_aadt','aadt_vn','aadt',clear_field_alias=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge major/minor roads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Merge minor roads.\n",
    "_ = arcpy.management.Merge(['centerlines_aadt',minor_roads],'centerlines_merged')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rename/Delete Fields\n",
    "For consistancy, readabiliy of the Python code, let's rename some of the fields and delete ones we don't care about"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['OBJECTID_1', 'Shape', 'Join_Count', 'TARGET_FID', 'OBJECTID', 'ADDR_SYS', 'CARTOCODE', 'FULLNAME', 'L_F_ADD', 'L_T_ADD', 'R_F_ADD', 'R_T_ADD', 'PREDIR', 'STREETNAME', 'STREETTYPE', 'SUFDIR', 'ALIAS1', 'ALIAS1TYPE', 'ALIAS2', 'ALIAS2TYPE', 'ACSALIAS', 'ACSNAME', 'ACSSUF', 'ADDR_QUAD', 'USPS_PLACE', 'ZIPLEFT', 'ZIPRIGHT', 'COFIPS', 'ONEWAY', 'SPEED', 'VERTLEVEL', 'CLASS', 'HWYNAME', 'DOT_RTNAME', 'DOT_RTPART', 'DOT_F_MILE', 'DOT_T_MILE', 'MODIFYDATE', 'GLOBALID', 'COLLDATE', 'ACCURACY', 'SOURCE', 'NOTES', 'COUNIQUE', 'SURFTYPE', 'SURFWIDTH', 'DSTRBWIDTH', 'LOCALFUNC', 'MAINTJURIS', 'FED_RDID', 'STATUS', 'ACCESS', 'USAGENOTES', 'DOT_RTID', 'DOT_FUNC', 'DOT_COFUND', 'LOCALID', 'aadt', 'Shape_Length']\n"
     ]
    }
   ],
   "source": [
    "print([f.name for f in arcpy.ListFields('centerlines_merged')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fields we don't need in the centerlines dataset for this\n",
    "delete_fields = [\n",
    "    'ADDR_SYS',\n",
    "    'CARTOCODE',\n",
    "    'FULLNAME',\n",
    "    'L_F_ADD',\n",
    "    'L_T_ADD',\n",
    "    'R_F_ADD',\n",
    "    'R_T_ADD',\n",
    "    'STREETNAME',\n",
    "    'ALIAS1',\n",
    "    'ALIAS1TYPE',\n",
    "    'ALIAS2',\n",
    "    'ALIAS2TYPE',\n",
    "    'ACSALIAS',\n",
    "    'ACSNAME',\n",
    "    'ACSSUF',\n",
    "    'ADDR_QUAD',\n",
    "    'USPS_PLACE', \n",
    "    'ZIPLEFT', \n",
    "    'ZIPRIGHT', \n",
    "    'COFIPS',\n",
    "    'VERTLEVEL', \n",
    "    'CLASS', \n",
    "    'HWYNAME', \n",
    "    'DOT_RTNAME',\n",
    "    'DOT_RTPART',\n",
    "    'DOT_F_MILE', \n",
    "    'DOT_T_MILE', \n",
    "    'MODIFYDATE',\n",
    "    'GLOBALID', \n",
    "    'COLLDATE', \n",
    "    'ACCURACY', \n",
    "    'SOURCE', \n",
    "    'NOTES', \n",
    "    'COUNIQUE', \n",
    "    'DSTRBWIDTH', \n",
    "    'LOCALFUNC',\n",
    "    'MAINTJURIS', \n",
    "    'FED_RDID', \n",
    "    'STATUS', \n",
    "    'ACCESS', \n",
    "    'USAGENOTES',\n",
    "    'DOT_RTID',\n",
    "    'DOT_FUNC', \n",
    "    'DOT_COFUND',\n",
    "    'LOCALID'\n",
    "]\n",
    "field_remap = {\n",
    "    'PREDIR':'pre_dir',\n",
    "    'STREETTYPE':'street_type',\n",
    "    'SUFDIR':'suf_dir',\n",
    "    'ONEWAY':'one_way',\n",
    "    'SPEED':'speed_limit',\n",
    "    'SURFTYPE':'surface_type',\n",
    "    'SURFWIDTH':'surface_width'\n",
    "}\n",
    "_ = arcpy.management.DeleteField('centerlines_merged',delete_fields)\n",
    "for field_name,remap_name in field_remap.items():\n",
    "    arcpy.management.AlterField('centerlines_merged',field_name,remap_name,clear_field_alias=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['OBJECTID_1', 'Shape', 'Join_Count', 'TARGET_FID', 'OBJECTID', 'pre_dir', 'street_type', 'suf_dir', 'one_way', 'speed_limit', 'surface_type', 'surface_width', 'aadt', 'Shape_Length']\n"
     ]
    }
   ],
   "source": [
    "print([f.name for f in arcpy.ListFields('centerlines_merged')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will reference *segment_id* as the unique identifier for a road segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Add a segment_id field to use instead of OBJECTID_1\n",
    "_ = arcpy.management.AddField('centerlines_merged','segment_id','Long')\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','segment_id','!OBJECTID_1!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Road Segment Spatial Features\n",
    "There are several fields to add to the data to enrich. Some will be calculated off of the geometries, some off of proximity to features in other datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Result './utah.gdb\\\\centerlines_merged'>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now we add some calculated fields:\n",
    "\n",
    "fields = [\n",
    "    ['sinuosity','Double'],\n",
    "    ['euclidean_length','Double'],\n",
    "    ['segment_length','Double'],\n",
    "    ['at_intersection','Short'],\n",
    "    ['near_billboard','Short'],\n",
    "    ['road_orient_approx','Double'],\n",
    "    ['near_major_road','Short']\n",
    "]\n",
    "\n",
    "_ = arcpy.management.AddFields('centerlines_merged',fields)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc Sinuosity\n",
    "code_block = \\\n",
    "'''\n",
    "import math\n",
    "def getSinuosity(shp):\n",
    "    x0 = shp.firstPoint.x\n",
    "    y0 = shp.firstPoint.y\n",
    "    x1 = shp.lastPoint.x\n",
    "    y1 = shp.lastPoint.y\n",
    "\n",
    "    euclid = math.sqrt((x0-x1)**2 + (y0-y1)**2)\n",
    "    length = shp.length\n",
    "    if euclid > 0:\n",
    "        return length/euclid\n",
    "    return 1.0\n",
    "'''\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','sinuosity','getSinuosity(!Shape!)',code_block=code_block)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc Euclidean Length\n",
    "code_block = \\\n",
    "'''\n",
    "import math\n",
    "def getEuclideanLength(shp):\n",
    "    x0 = shp.firstPoint.x\n",
    "    y0 = shp.firstPoint.y\n",
    "    x1 = shp.lastPoint.x\n",
    "    y1 = shp.lastPoint.y\n",
    "\n",
    "    euclid = math.sqrt((x0-x1)**2 + (y0-y1)**2)\n",
    "    return euclid\n",
    "'''\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','euclidean_length','getEuclideanLength(!Shape!)',code_block=code_block)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc segment_length\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','segment_length','!Shape_Length!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc Road Orientation\n",
    "code_block = \\\n",
    "'''\n",
    "import math\n",
    "def getOrientation(shp):\n",
    "    x0 = shp.firstPoint.x\n",
    "    y0 = shp.firstPoint.y\n",
    "    x1 = shp.lastPoint.x\n",
    "    y1 = shp.lastPoint.y\n",
    "\n",
    "    angle = math.atan2(y1-y0,x1-x0)\n",
    "    if angle < 0:\n",
    "        angle += math.pi\n",
    "    return angle\n",
    "'''\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','road_orient_approx','getOrientation(!Shape!)',code_block=code_block)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc at_intersection --Requires ArcGIS Pro Advanced License as it uses 'Near' tool \n",
    "_ = arcpy.analysis.Near('centerlines_merged',\n",
    "                        intersections_url,\n",
    "                        search_radius='150 feet')\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','at_intersection','!NEAR_DIST! > 0')\n",
    "_ = arcpy.management.DeleteField('centerlines_merged',['NEAR_DIST','NEAR_FID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc near_billboard\n",
    "_ = arcpy.analysis.Near('centerlines_merged',\n",
    "                        billboards_url,\n",
    "                        search_radius='300 feet')\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','near_billboard','!NEAR_DIST! > 0')\n",
    "_ = arcpy.management.DeleteField('centerlines_merged',['NEAR_DIST','NEAR_FID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calc near_major_road\n",
    "_ = arcpy.analysis.Near('centerlines_merged',\n",
    "                        'centerlines',\n",
    "                        search_radius='50 feet')\n",
    "_ = arcpy.management.CalculateField('centerlines_merged','near_major_road','!NEAR_DIST! > 0')\n",
    "_ = arcpy.management.DeleteField('centerlines_merged',['NEAR_DIST','NEAR_FID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load collisions\n",
    "\n",
    "# Projected Coords UTM NAD1983 Zone 12N\n",
    "in_sr = arcpy.SpatialReference(26912)\n",
    "\n",
    "_ = arcpy.management.MakeXYEventLayer(collisions_csv,\n",
    "        out_layer='collisions',\n",
    "        in_x_field='Longitude',\n",
    "        in_y_field='Latitude',\n",
    "        spatial_reference=in_sr)\n",
    "\n",
    "# Projected Coords Web Mercator\n",
    "out_sr = arcpy.Describe('centerlines_merged').spatialReference\n",
    "\n",
    "_ = arcpy.management.Project('collisions','collisions_projected',out_sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assign Weather Stations to Segments\n",
    "We need to create the Thiessen polygons for the weather station"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "in_sr = arcpy.SpatialReference(4326)\n",
    "_ = arcpy.management.MakeXYEventLayer(weather_stations_csv,\n",
    "        out_layer='weather_station_eventlayer',\n",
    "        in_x_field='LON',\n",
    "        in_y_field='LAT',\n",
    "        spatial_reference=in_sr)\n",
    "res = arcpy.management.CopyFeatures('weather_station_eventlayer','weather_station_locations',)\n",
    "_ = arcpy.analysis.CreateThiessenPolygons(res,'weather_station_polygons',fields_to_copy='ALL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge weather station polygons with the roads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "field_mappings = arcpy.FieldMappings()\n",
    "field_mappings.addTable('centerlines_merged')\n",
    "features_mapping = arcpy.FieldMap()\n",
    "features_mapping.addInputField('weather_station_polygons','station_id')\n",
    "field_mappings.addFieldMap(features_mapping)\n",
    "\n",
    "_ = arcpy.analysis.SpatialJoin('centerlines_merged',\n",
    "                               'weather_station_polygons',\n",
    "                               'centerlines_merged_with_weather',\n",
    "                               join_type='KEEP_ALL',\n",
    "                               match_option='closest',\n",
    "                               search_radius='100 miles',\n",
    "                               field_mapping=field_mappings)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Result './utah.gdb\\\\closest_weather_stations'>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_ = arcpy.analysis.GenerateNearTable('centerlines_merged',\n",
    "                                     'weather_station_locations',\n",
    "                                     'closest_weather_stations',\n",
    "                                     closest=False,\n",
    "                                     closest_count=4)\n",
    "arcpy.management.AlterField(_,'IN_FID','segment_id',clear_field_alias=True)\n",
    "arcpy.management.AlterField(_,'NEAR_FID','station_id',clear_field_alias=True)\n",
    "arcpy.management.AlterField(_,'NEAR_DIST','distance',clear_field_alias=True)\n",
    "arcpy.management.AlterField(_,'NEAR_RANK','rank',clear_field_alias=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_ = arcpy.management.CopyFeatures('centerlines_merged_with_weather','static_features')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatially Assign Collisions to Roads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "field_mappings = arcpy.FieldMappings()\n",
    "field_mappings.addTable('collisions_projected')\n",
    "\n",
    "features_mapping = arcpy.FieldMap()\n",
    "features_mapping.addInputField('static_features','segment_id')\n",
    "field_mappings.addFieldMap(features_mapping)\n",
    "\n",
    "_ = arcpy.analysis.SpatialJoin('collisions_projected',\n",
    "                               'static_features',\n",
    "                               'collisions_joined',\n",
    "                               join_type='KEEP_COMMON', # Right join\n",
    "                               match_option='closest',\n",
    "                               search_radius='50 feet',\n",
    "                               field_mapping=field_mappings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}