417 lines
12 KiB
Plaintext
417 lines
12 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from matplotlib.patches import Polygon\n",
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"import fiona\n",
|
|
"from shapely.geometry import shape\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import os\n",
|
|
"from collections import defaultdict, OrderedDict \n",
|
|
"\n",
|
|
"## Fiona, ogr2ogr, pyplot demo\n",
|
|
"## Live 8.5 * darkblue-b\n",
|
|
"##"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"`Demo Area Analysis in IPython Notebook`\n",
|
|
"\n",
|
|
" * Use North Carolina sample data supplied on the OSGeo Live \n",
|
|
" * Extract an area of interest\n",
|
|
" * Extract a subset of a state soils record based on a bounding box (BBOX)\n",
|
|
" * Plot geometry graphically\n",
|
|
" * Buffer the geometry by a fixed amount; calculate the areas of intersection for soil types\n",
|
|
" * Produce a chart of the totals with pyplot"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"infile = '/usr/local/share/data/north_carolina/shape/urbanarea.shp'\n",
|
|
"infile_soils = '/usr/local/share/data/north_carolina/shape/soils_general.shp'\n",
|
|
"\n",
|
|
"out_farmville = '/tmp/farmville_shp'\n",
|
|
"out_farmville_soils = '/tmp/farmville_soil_shp'"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## -- Call the gdal/ogr tool suite to extract and subset the original statewide dataset\n",
|
|
"## use a feature of IPython Notebook to pass python variables to a command line\n",
|
|
"\n",
|
|
"!rm -rf $out_farmville\n",
|
|
"!ogr2ogr -f 'ESRI Shapefile' $out_farmville $infile -where \"NAME='Farmville'\" -dim 2 -a_srs EPSG:3358\n",
|
|
"\n",
|
|
"##-- advanced example -- raleigh is multiple polygons\n",
|
|
"#!rm -rf /tmp/raleigh1_shp\n",
|
|
"#!ogr2ogr -f 'ESRI Shapefile' /tmp/raleigh1_shp $infile -where \"NAME='Raleigh'\" -dim 2 -a_srs EPSG:3358\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## -- use the IPython notebook to get help on py module Fiona\n",
|
|
"# fiona?\n",
|
|
"# fiona.open?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## open and inspect a Shapefile using Fiona\n",
|
|
"## a POLYGON record could be very long, so dont print the record w/ geometry\n",
|
|
"with fiona.open( out_farmville ) as f:\n",
|
|
" crs = f.crs\n",
|
|
" print( 'CRS:',crs)\n",
|
|
" rec = f.next()\n",
|
|
" #print rec\n",
|
|
" print( 'SHAPELY REC:',rec.keys() )\n",
|
|
" print( 'SHAPEFILE FLDS: ',rec['properties'].keys() )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"with fiona.open( out_farmville,'r') as inp:\n",
|
|
" #print dir(inp) ## whats inside this object?\n",
|
|
" #print inp.bounds\n",
|
|
" \n",
|
|
" ## take the bounding box of the area of interest\n",
|
|
" ## add 300 meters on each side for aesthetics\n",
|
|
" ## (left, bottom, right, top)\n",
|
|
" left_bnds = inp.bounds[0] - 300\n",
|
|
" bottom_bnds = inp.bounds[1] - 300\n",
|
|
" right_bnds = inp.bounds[2] + 300\n",
|
|
" top_bnds = inp.bounds[3] + 300\n",
|
|
"\n",
|
|
"## echo one variable to sanity check\n",
|
|
"#left_bnds"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## ogr2ogr ... -clipsrc [xmin ymin xmax ymax] ...\n",
|
|
"\n",
|
|
"!rm -rf $out_farmville_soils\n",
|
|
"!ogr2ogr -f 'ESRI Shapefile' $out_farmville_soils $infile_soils -clipsrc $left_bnds $bottom_bnds $right_bnds $top_bnds -dim 2 -a_srs EPSG:3358\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## ----------------------------------------\n",
|
|
"## plot_multipolygon function\n",
|
|
"## Author: Kelsey Jordahl, Enthought\n",
|
|
"## Scipy 2013 geospatial tutorial\n",
|
|
"\n",
|
|
"def plot_polygon(ax, poly, color='black'):\n",
|
|
" a = np.asarray(poly.exterior)\n",
|
|
" ax.add_patch(Polygon(a, facecolor=color, alpha=0.5))\n",
|
|
" ax.plot(a[:, 0], a[:, 1], color=color)\n",
|
|
"\n",
|
|
"def plot_multipolygon(ax, geom, color='red'):\n",
|
|
" \"\"\" Can safely call with either Polygon or Multipolygon geometry\n",
|
|
" \"\"\"\n",
|
|
" if geom.type == 'Polygon':\n",
|
|
" plot_polygon(ax, geom, color)\n",
|
|
" elif geom.type == 'MultiPolygon':\n",
|
|
" for poly in geom.geoms:\n",
|
|
" plot_polygon(ax, poly, color)\n",
|
|
"## ----------------------------------------"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"soil_colors = ['green', 'brown', '#ceFFde']\n",
|
|
"\n",
|
|
"dst_geoms = OrderedDict()\n",
|
|
"src_geoms = OrderedDict()\n",
|
|
"\n",
|
|
"cnt = 0\n",
|
|
"with fiona.open( out_farmville ) as f:\n",
|
|
" for rec in f:\n",
|
|
" src_geoms[cnt] = shape(rec['geometry'])\n",
|
|
" cnt += 1\n",
|
|
"\n",
|
|
"cnt = 0\n",
|
|
"with fiona.open( out_farmville_soils ) as f:\n",
|
|
" for rec in f:\n",
|
|
" ## check the geometry type if desired\n",
|
|
" ## intersections may result in POINT or LINESTRING\n",
|
|
" if rec['geometry']['type'] != 'Polygon':\n",
|
|
" continue\n",
|
|
" gsl = rec['properties']['GSL_NAME']\n",
|
|
" dst_geoms[gsl] = shape(rec['geometry'])\n",
|
|
" cnt += 1\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\n",
|
|
"fig, ax = plt.subplots(figsize=(11,11))\n",
|
|
"plt.title(\"Soil Class and Urban Area via N. Carolina Sample Dataset (CC-by-SA) OSGeo\")\n",
|
|
"\n",
|
|
"cnt = 0\n",
|
|
"for key in dst_geoms :\n",
|
|
" ## cnt mod (number of colors) is always a safe index into colors[]\n",
|
|
" color = soil_colors[ cnt%len(soil_colors) ]\n",
|
|
" plot_multipolygon(ax, dst_geoms[key], color=color)\n",
|
|
" cnt += 1\n",
|
|
" \n",
|
|
"cnt = 0\n",
|
|
"color = 'gray'\n",
|
|
"with fiona.open( out_farmville ) as f:\n",
|
|
" for rec in f:\n",
|
|
" plot_multipolygon(ax, src_geoms[cnt], color=color)\n",
|
|
" cnt += 1\n",
|
|
"\n",
|
|
"#ax.add_patch(Polygon( src_geoms[0].centroid, facecolor='black', alpha=0.5))\n",
|
|
"\n",
|
|
"labels = ax.get_xticklabels() \n",
|
|
"for label in labels: \n",
|
|
" label.set_rotation(90) \n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## show all the variables defined so far\n",
|
|
"%whos"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## what geomtries have been created ?\n",
|
|
"print( 'src_geoms len: ',len(src_geoms) )\n",
|
|
"print( 'dst_geoms len: ',len(dst_geoms ) )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## buffer the urbanarea by N meters\n",
|
|
"## save the result in the source-geometry list\n",
|
|
"\n",
|
|
"src_geom_buffered = src_geoms[0].buffer(166)\n",
|
|
"src_geoms[1] = src_geom_buffered\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, ax = plt.subplots(figsize=(11,11))\n",
|
|
"\n",
|
|
"cnt = 0\n",
|
|
"for key in dst_geoms :\n",
|
|
" color = soil_colors[ cnt%len(soil_colors) ]\n",
|
|
" plot_multipolygon(ax, dst_geoms[key], color=color)\n",
|
|
" cnt += 1\n",
|
|
"\n",
|
|
"plot_multipolygon(ax, src_geoms[1], color='gray')\n",
|
|
"plot_multipolygon(ax, src_geoms[0], color='gray')\n",
|
|
" \n",
|
|
"labels = ax.get_xticklabels() \n",
|
|
"for label in labels: \n",
|
|
" label.set_rotation(90) \n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## take the intersection of each soil area against the buffered poly of interest\n",
|
|
"## store into a third convenient list\n",
|
|
"res_geoms = OrderedDict()\n",
|
|
"for key in dst_geoms:\n",
|
|
" res_geoms[key] = src_geom_buffered.intersection( dst_geoms[key] )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, ax = plt.subplots(figsize=(11,11))\n",
|
|
"\n",
|
|
"cnt = 0\n",
|
|
"for key in res_geoms :\n",
|
|
" color = soil_colors[ cnt%len(soil_colors) ]\n",
|
|
" plot_multipolygon(ax, res_geoms[key], color=color)\n",
|
|
" cnt += 1\n",
|
|
" \n",
|
|
"labels = ax.get_xticklabels() \n",
|
|
"for label in labels: \n",
|
|
" label.set_rotation(90) \n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## define a lookup table of mock attributes of each soil type\n",
|
|
"## (substitute a real data source here)\n",
|
|
"\n",
|
|
"soils_lkup = {\n",
|
|
" 'NC038' : { 'TEXTURE_CLAY':0.53, 'TEXTURE_SILT':0.33, 'TEXTURE_SAND':0.44 },\n",
|
|
" 'NC035' : { 'TEXTURE_CLAY':0.70, 'TEXTURE_SILT':0.33, 'TEXTURE_SAND':0.44 },\n",
|
|
" 'NC034' : { 'TEXTURE_CLAY':0.23, 'TEXTURE_SILT':0.74, 'TEXTURE_SAND':0.44 }\n",
|
|
"}\n",
|
|
"\n",
|
|
"## get a rough total area for display purposes\n",
|
|
"sum_areas = 0.0\n",
|
|
"for key in res_geoms:\n",
|
|
" sum_areas += int(res_geoms[key].area)\n",
|
|
"\n",
|
|
"## record the area and percentage area in a convenient dictionary \n",
|
|
"tdd = {}\n",
|
|
"for key in res_geoms:\n",
|
|
" tdd[key] = soils_lkup[key]\n",
|
|
" tdd[key]['area'] = int(res_geoms[key].area)\n",
|
|
" tdd[key]['perc'] = int((res_geoms[key].area/sum_areas)*100)/100.0\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## Now all attributes for visualization are available in a single dict\n",
|
|
"## in a larger system this could be delivered for serving graphical reports\n",
|
|
"tdd"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"## What matplotlib pylot graphs are available ?\n",
|
|
"## http://matplotlib.org/api/pyplot_api.html\n",
|
|
"##\n",
|
|
"## Use IPython built-in help to discover pie chart attributes\n",
|
|
"\n",
|
|
"plt.pie?\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# The slices will be ordered and plotted counter-clockwise.\n",
|
|
"labels = [ 'NC034', 'NC035', 'NC038' ]\n",
|
|
"sizes = [ tdd['NC034']['perc']*100, tdd['NC035']['perc']*100, tdd['NC038']['perc']*100 ]\n",
|
|
"pie_colors = ['lightcoral', '#eeefee', 'yellowgreen']\n",
|
|
"# only \"explode\" the 2nd slice (i.e. 'NC035')\n",
|
|
"explode = (0, 0.1, 0)\n",
|
|
"\n",
|
|
"plt.pie(sizes, explode=explode, labels=labels, colors=pie_colors,\n",
|
|
" autopct='%1.1f%%', shadow=True, startangle=90)\n",
|
|
"# Set aspect ratio to be equal so that pie is drawn as a circle.\n",
|
|
"plt.axis('equal')\n",
|
|
"\n",
|
|
"plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"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.8.5"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 1
|
|
}
|