Skip to content

Instantly share code, notes, and snippets.

@sjsrey

sjsrey/airports.ipynb

Last active August 29, 2015 14:23
Show Gist options
  • Save sjsrey/8f68a4518521b95227ce to your computer and use it in GitHub Desktop.
Save sjsrey/8f68a4518521b95227ce to your computer and use it in GitHub Desktop.
Download ZIP
Exploring PySAL arc based KDTrees for nearest neighbor weights
Display the source blob
Display the rendered blob
Raw
airports.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exploring PySAL's KDTree for nearest neighbor weights on the sphere\n",
"\n",
"Files used in this notebook are available at https://gist.github.com/sjsrey/f5212448c4053e27e329\n",
"\n",
"Author: Serge Rey <sjsrey@gmail.com>"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import pysal as ps\n",
"import numpy as np\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#files available at https://gist.github.com/sjsrey/f5212448c4053e27e329\n",
"shpf = ps.open(\"airports.shp\")\n",
"points = []\n",
"for point in shpf:\n",
" points.append(point)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"6977"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(points)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"6977"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"point.id"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(-4.1749, 55.5141)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"point"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"npoints = np.array(points)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 145.391881, -6.081689],\n",
" [ 145.7887 , -5.207083],\n",
" [ 144.295861, -5.826789],\n",
" ..., \n",
" [-121.1825 , 51.3833 ],\n",
" [-121.1958 , 51.4412 ],\n",
" [ -4.1749 , 55.5141 ]])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"npoints"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"tree = ps.cg.KDTree(npoints, distance_metric=\"Arc\",radius=6371.0)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nn3 = tree.query(tree.data,k=4)#self is always closest so bump k by 1"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([[ 0. , 46.97452698, 106.70511125, 109.83659115],\n",
" [ 0. , 106.70511125, 128.23847708, 163.18754307],\n",
" [ 0. , 16.22095444, 48.81862711, 77.8333434 ],\n",
" ..., \n",
" [ 0. , 6.50392731, 91.63902287, 107.25361766],\n",
" [ 0. , 6.50392731, 97.50720006, 101.42195232],\n",
" [ 0. , 25.93108079, 39.0012831 , 39.32575656]]),\n",
" array([[ 0, 4901, 1, 6368],\n",
" [ 1, 0, 4901, 6368],\n",
" [ 2, 6368, 4884, 4901],\n",
" ..., \n",
" [6974, 6975, 77, 160],\n",
" [6975, 6974, 77, 160],\n",
" [6976, 532, 5962, 6240]]))"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"['airport_id',\n",
" 'name',\n",
" 'city',\n",
" 'country',\n",
" 'iata_faa',\n",
" 'iaco',\n",
" 'latitude',\n",
" 'longitude',\n",
" 'altitude',\n",
" 'zone',\n",
" 'dst']"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dbf = ps.open(\"airports.dbf\")\n",
"dbf.header"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"airport_id = dbf.by_col('airport_id')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1992"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id.index(2033)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0. , 1349.0885623 , 1351.26578069, 2446.28820941])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3[0][1992]"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"name = dbf.by_col('name')\n",
"country = dbf.by_col('country')"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"('South Pole Station', 'Antarctica')"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"name[1992], country[1992]"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([1992, 1997, 6520, 5308])"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3[1][1992]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"7947"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id[6520]"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2038"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id[1997]"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(6437, 'Mcmurdo Station')"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id[5308], name[1997]"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"country = dbf.by_col('country')\n",
"name = dbf.by_col('name')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Antarctica Mcmurdo Station\n",
"Antarctica Ice Runway-McMurdo Sstation\n",
"Antarctica South Shetland\n"
]
}
],
"source": [
"for i in nn3[1][1992][1:]:\n",
" print country[i], name[i]"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4971"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id.index(6098)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Russia Ugolny Airport\n",
"Russia Provideniya Bay\n",
"United States Gambell Airport\n",
"United States Savoonga Airport\n"
]
}
],
"source": [
"for i in nn3[1][4971]:\n",
" print country[i], name[i]"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([4971, 2852, 5517, 5508])"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3[1][4971]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"6715"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"airport_id[5517]"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0. , 432.12388349, 519.07360721, 580.09537777])"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3[0][4971]"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# weights based on arc kdtree"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"w3 = ps.knnW(tree,k=3)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[2852, 5517, 5508]"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w3.neighbors[4971]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"'Russia'"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"country[4971]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[('Russia', 'Provideniya Bay'),\n",
" ('United States', 'Gambell Airport'),\n",
" ('United States', 'Savoonga Airport')]"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[(country[j],name[j]) for j in w3.neighbors[4971]]"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"['United States', 'Russia', 'United States']"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[country[j] for j in w3.neighbors[5517]]"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"nn3d = nn3[0]"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(6977, 4)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3d.shape"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0. , 2599.29658612, 2907.40848333, 2925.68682231])"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nn3d.max(axis=0)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"maxnnd = nn3d.max(axis=0)[1]"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([2587]),)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.where(nn3d[:,1]==maxnnd)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"airport_id[2587]\n",
"city = dbf.by_col('city')"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"('Chile', 'Mataveri Intl', 'Easter Island')"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# most isolated airport from nearest neighbor\n",
"country[2587], name[2587], city[2587] \n"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# probably why the local population went extinct ;->"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment