{ "cells": [ { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "a20290f5-2a26-4622-9e8b-ce94cc1dcb1a" }, "slideshow": { "slide_type": "slide" } }, "source": [ "# Python & GeoComputation\n", "Longzhu Shen Spatial Ecology Jun 2019\n", "\n", "Code avalability at\n", "\n", " wget https://github.com/selvaje/spatial-ecology-codes/blob/master/docs/source/PYTHON/02_Geo_Python.ipynb" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Configuration\n", "\n", "Adjust the 'ast_note_interactivity' setting to print all versus last expression in Jupyter cell." ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [], "source": [ "#from IPython.core.interactiveshell import InteractiveShell\n", "#InteractiveShell.ast_node_interactivity = \"last_expr\" #\"all\" " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Numpy Intro\n", "\n", "Numpy is a foundational package for data analysis in Python. It enables extremely efficient mathematical operations on large data sets. It also can be used to store and manipulate tabular data." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Creating N-dimensional arrays using NumPy\n", "\n", "- There are many ways to create N-dimensional arrays." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Create 2X3 array initialized to all zeroes & check the dimensions of the array" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0.],\n", " [0., 0., 0.]])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.zeros((2,3))\n", "a" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 3)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.shape" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Create array using arange function" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 2, 3, 4, 5])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.arange(6)\n", "a" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Create array initialized by list of lists" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2],\n", " [3, 4, 5]])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([[0,1,2],[3,4,5]])\n", "a" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Get values from N-dimensional array" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0,0] " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "slicing" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 2])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0,:] " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([2, 5])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[:,-1]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[1, 2],\n", " [4, 5]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0:2,1:3]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Modifying N-dimensional arrays" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2],\n", " [3, 4, 5]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[25, 1, 2],\n", " [ 3, 4, 5]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0,0] = 25.0\n", "a" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[10, 11, 12],\n", " [ 3, 4, 5]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0,:] = np.array([10,11,12])\n", "a" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[10, 11, 20],\n", " [ 3, 4, 21]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[:,-1] = np.array([20,21])\n", "a" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[10, 10, 11],\n", " [ 3, 20, 21]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0:2,1:3] = np.array([[10,11],[20,21]])\n", "a" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Broadcasting\n", "\n", "\n", "- \"Describes how numpy treats arrays with different shapes during arithmetic operations\"\n", "- Subject to constraints, smaller array is 'broadcast' across larger array so they have compatible shapes\n", "- Looping occurs in C instead of Python - fast!\n", "\n", "More info [here](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[10, 10, 10],\n", " [ 3, 4, 21]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0,:] = 10.0\n", "a" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[30, 31, 32],\n", " [30, 31, 32]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[:,:] = np.array([30,31,32])\n", "a" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[40],\n", " [41]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tmp = np.array([40,41]).reshape(2,1)\n", "tmp" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[40, 40, 40],\n", " [41, 41, 41]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[:,:] = tmp\n", "a" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[ 40, 100, 100],\n", " [ 41, 100, 100]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a[0:2,1:3] = 100.0\n", "a" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Math Operations on Arrays\n", "\n", "- Arithmetic of arrays vs. scalars applies scalar to all elements of array\n", "- Arithmetic of array vs. array operates elementwise" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 2, 3])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.arange(4)\n", "a" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.5, 1.5, 2.5, 3.5])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a + 0.5" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 2, 4, 6])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a + a" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 4, 9])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a * a" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "14" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dot(a, a)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([[0, 0, 0, 0],\n", " [0, 1, 2, 3],\n", " [0, 2, 4, 6],\n", " [0, 3, 6, 9]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dot(a.reshape(4,1), a.reshape(1,4))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Working with Tabular Data" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "solar_system_file = 'solar_system_abbr.csv' \n", "solar_system_data = np.genfromtxt(solar_system_file, delimiter=',', skip_header=1,\n", " dtype=['S10', 'S20', 'i4', 'f4' , 'f4'])" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "array([(b'Sun', b'Star', 0, 1.3920e+06, 3.33e+05),\n", " (b'Mercury', b'Terrestrial planet', 1, 4.8780e+03, 5.50e-02),\n", " (b'Venus', b'Terrestrial planet', 2, 1.2104e+04, 8.15e-01),\n", " (b'Earth', b'Terrestrial planet', 3, 1.2756e+04, 1.00e+00),\n", " (b'Mars', b'Terrestrial planet', 4, 6.7870e+03, 1.07e-01),\n", " (b'Jupiter', b'Gas giant', 6, 1.4280e+05, 3.18e+02),\n", " (b'Saturn', b'Gas giant', 7, 1.2000e+05, 9.50e+01),\n", " (b'Uranus', b'Gas giant', 8, 5.1118e+04, 1.50e+01),\n", " (b'Neptune', b'Gas giant', 9, 4.9528e+04, 1.70e+01),\n", " (b'Ceres', b'Dwarf planet', 5, 9.7460e+02, 1.60e-04),\n", " (b'Pluto', b'Dwarf planet', 10, 2.3000e+03, 2.00e-03),\n", " (b'Haumea', b'Dwarf planet', 11, 1.3000e+03, 7.00e-04),\n", " (b'Makemake', b'Dwarf planet', 12, 1.4200e+03, 6.70e-04),\n", " (b'Eris', b'Dwarf planet', 13, 2.3260e+03, 2.80e-03)],\n", " dtype=[('f0', 'S10'), ('f1', 'S20'), ('f2', ')" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Retrieve a file \n", "import urllib.request\n", "url = \"http://spatial-ecology.net/dokuwiki/lib/exe/fetch.php?media=wiki:python:hancock.zip\"\n", "fileName = \"hancock.zip\"\n", "urllib.request.urlretrieve(url, fileName)\n", " " ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b'time,latitude,longitude,depth,mag,magType,nst,gap,dmin,rms,net,id,updated,place,type,horizontalError,depthError,magError,magNst,status,locationSource,magSource\\n'\n", "b'2019-05-31T17:17:57.190Z,34.0203323,-117.6094971,19.12,1.96,ml,13,134,0.1116,0.28,ci,ci38613536,2019-05-31T17:20:00.408Z,\"4km SSE of Ontario, CA\",earthquake,2.14,1.35,0.412,2,automatic,ci,ci\\n'\n", "b'2019-05-31T17:05:15.100Z,33.4946667,-116.795,1.59,0.49,ml,19,57,0.01821,0.23,ci,ci38613504,2019-05-31T17:08:56.963Z,\"9km NE of Aguanga, CA\",earthquake,0.45,0.67,0.127,10,automatic,ci,ci\\n'\n", "b'2019-05-31T16:58:50.720Z,33.4896667,-116.7948333,2.37,1.59,ml,56,36,0.02311,0.2,ci,ci38613480,2019-05-31T17:10:18.350Z,\"8km NE of Aguanga, CA\",earthquake,0.19,0.34,0.151,27,automatic,ci,ci\\n'\n", "b'2019-05-31T16:39:45.470Z,45.867,-111.361,4.94,1.57,ml,13,88,0.302,0.13,mb,mb80346349,2019-05-31T17:05:26.100Z,\"2km WNW of Manhattan, Montana\",earthquake,0.39,8.54,0.255,6,reviewed,mb,mb\\n'\n", "b'2019-05-31T16:37:53.690Z,33.5103333,-116.5108333,11.17,0.7,ml,25,82,0.04997,0.18,ci,ci38613456,2019-05-31T16:41:30.754Z,\"16km ESE of Anza, CA\",earthquake,0.3,0.53,0.161,16,automatic,ci,ci\\n'\n", "b'2019-05-31T16:36:34.830Z,34.0466667,-117.5038333,2.53,1.07,ml,31,64,0.06375,0.2,ci,ci38613448,2019-05-31T16:40:14.360Z,\"4km NNW of Glen Avon, CA\",earthquake,0.32,0.47,0.168,30,automatic,ci,ci\\n'\n", "b'2019-05-31T16:22:07.250Z,38.8214989,-122.8526688,2.31,1.13,md,20,64,0.006873,0.03,nc,nc73190601,2019-05-31T17:06:03.673Z,\"10km WNW of The Geysers, CA\",earthquake,0.23,0.56,0.03,4,automatic,nc,nc\\n'\n" ] } ], "source": [ "# Read data from a url - in this case recent USGS earthquakes\n", "import urllib.request\n", "url = \"http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_hour.csv\"\n", "earthquakes = urllib.request.urlopen(url)\n", "\n", "# Read the first two earthquakes\n", "print(earthquakes.readline())\n", "print(earthquakes.readline())\n", "\n", "# Iterate through the rest\n", "for record in earthquakes: \n", " print (record)" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "e675a434-eb22-4d12-a159-d6464ec9d916" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## Zipping Files" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Unzip a shapefile \n", "import zipfile\n", "zip = open(\"hancock.zip\", \"rb\") # rb = 'read binary'; more info @ https://docs.python.org/2/tutorial/inputoutput.html\n", "zipShape = zipfile.ZipFile(zip)\n", "for fileName in zipShape.namelist():\n", " out = open(fileName, \"wb\")\n", " out.write(zipShape.read(fileName))\n", " out.close()" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Add a shapefile to a tar archive\n", "import tarfile\n", "tar = tarfile.open(\"hancock.tar.gz\", \"w:gz\")\n", "tar.add(\"hancock.shp\")\n", "tar.add(\"hancock.shx\")\n", "tar.add(\"hancock.dbf\")\n", "tar.close()\n", "\n", "tar" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Extract a shapefile from a gzipped tar archive\n", "import tarfile\n", "tar = tarfile.open(\"hancock.tar.gz\", \"r:gz\")\n", "tar.extractall()\n", "tar.close()\n", " \n", "tar" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "e675a434-eb22-4d12-a159-d6464ec9d916" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## Raster Images" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### Loading Files" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3\n", "2592\n", "2693\n" ] } ], "source": [ "# Open a raster with gal\n", "import gdal\n", "\n", "raster = gdal.Open(\"SatImage.tif\")\n", "\n", "print (raster.RasterCount)\n", "print (raster.RasterXSize)\n", "print (raster.RasterYSize)\n", " " ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ " >" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Numpy/gdalnumeric - Read an image, extract a band, save a new image\n", "from osgeo import gdalnumeric\n", "\n", "srcArray = gdalnumeric.LoadFile(\"SatImage.tif\")\n", "band1 = srcArray[0]\n", "gdalnumeric.SaveArray(band1, \"band1.jpg\", format=\"JPEG\")\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### Rasterization" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Rasterize a shapefile with PIL\n", "import shapefile\n", "from PIL import Image, ImageDraw\n", "\n", "r = shapefile.Reader(\"hancock.shp\")\n", "xdist = r.bbox[2] - r.bbox[0]\n", "ydist = r.bbox[3] - r.bbox[1]\n", "iwidth = 400\n", "iheight = 600\n", "xratio = iwidth/xdist\n", "yratio = iheight/ydist\n", "pixels = []\n", "for x,y in r.shapes()[0].points:\n", " px = int(iwidth - ((r.bbox[2] - x) * xratio))\n", " py = int((r.bbox[3] - y) * yratio)\n", " pixels.append((px,py))\n", "img = Image.new(\"RGB\", (iwidth, iheight), \"white\")\n", "draw = ImageDraw.Draw(img)\n", "draw.polygon(pixels, outline=\"rgb(203, 196, 190)\", fill=\"rgb(198, 204, 189)\")\n", "img.save(\"hancock.png\")" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create a pdf with a map image\n", " \n", "import fpdf\n", " \n", "# PDF constructor:\n", "# Portrait, millimeter units, A4 page size\n", "pdf=fpdf.FPDF(\"P\", \"mm\", \"A4\")\n", "# create a new page\n", "pdf.add_page()\n", "# Set font: arial, bold, size 20\n", "pdf.set_font('Arial','B',20)\n", "# Create a new cell: 160 x 25mm, text contents, no border, centered\n", "pdf.cell(160,25,'Hancock County Boundary', border=0, align=\"C\")\n", "pdf.image(\"hancock.png\",25,50,150,160)\n", "pdf.output('map.pdf','F')\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### Raster Processing" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ " >" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"Swap bands in a raster satellite image\"\"\"\n", "\n", "# Module within the GDAL python package \n", "import gdalnumeric\n", "\n", "# name of our source image\n", "src = \"FalseColor.tif\"\n", "\n", "# load the source image into an array\n", "arr = gdalnumeric.LoadFile(src)\n", "\n", "# swap bands 1 and 2 for a natural color image.\n", "# We will use numpy \"advanced slicing\" to reorder the bands.\n", "# Using the source image\n", "gdalnumeric.SaveArray(arr[[1,0,2],:], \"swap.tif\", format=\"GTiff\", prototype=src)\n" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ " >" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"Perform a simple difference image change detection on a \n", "'before' and 'after' image.\"\"\"\n", "import gdal, gdalnumeric\n", "import numpy as np\n", " \n", "# \"Before\" image\n", "im1 = \"before.tif\"\n", "\n", "# \"After\" image\n", "im2 = \"after.tif\"\n", "\n", "# Load before and after into arrays\n", "ar1 = gdalnumeric.LoadFile(im1).astype(np.int8) \n", "ar2 = gdalnumeric.LoadFile(im2)[1].astype(np.int8)\n", "\n", "# Perform a simple array difference on the images\n", "diff = ar2 - ar1\n", "\n", "# Set up our classification scheme to try\n", "# and isolate significant changes\n", "classes = np.histogram(diff, bins=5)[1]\n", "\n", "# The color black is repeated to mask insignificant changes\n", "lut = [[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,255,0],[255,0,0]]\n", "\n", "# Starting value for classification\n", "start = 1\n", "\n", "# Set up the output image\n", "rgb = np.zeros((3, diff.shape[0], diff.shape[1],), np.int8) \n", "\n", "# Process all classes and assign colors\n", "for i in range(len(classes)):\n", " mask = np.logical_and(start <= diff, diff <= classes[i])\n", " for j in range(len(lut[i])):\n", " rgb[j] = np.choose(mask, (rgb[j], lut[i][j]))\n", " start = classes[i]+1\n", "\n", " # Save the output image\n", "gdalnumeric.SaveArray(rgb, \"change.tif\", format=\"GTiff\", prototype=im2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "e675a434-eb22-4d12-a159-d6464ec9d916" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## Shapefiles" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 1.0 First\n", "3.0 1.0 Second\n", "4.0 3.0 Third\n", "2.0 2.0 Fourth\n", "0.0 0.0 Appended\n" ] } ], "source": [ "#Examine a shapefile with ogr \n", " \n", "from osgeo import ogr\n", "\n", "# open the shapefile\n", "shp = ogr.Open(\"point.shp\")\n", "\n", "# Get the layer\n", "layer = shp.GetLayer()\n", "\n", "# Loop through the features and print information about them\n", "for feature in layer:\n", " geometry = feature.GetGeometryRef()\n", " print (geometry.GetX(), geometry.GetY(), feature.GetField(\"FIRST_FLD\"))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 1.0 First\n", "3.0 1.0 Second\n", "4.0 3.0 Third\n", "2.0 2.0 Fourth\n", "0.0 0.0 Appended\n" ] } ], "source": [ "# Examine a shapefile with pyshp \n", "import shapefile\n", " \n", "shp = shapefile.Reader(\"point\")\n", "for feature in shp.shapeRecords():\n", " point = feature.shape.points[0]\n", " rec = feature.record[0]\n", " print (point[0], point[1], rec)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "\"\"\"Choropleth thematic map\"\"\"\n", "import math\n", "import shapefile\n", "from PIL import Image, ImageDraw\n", "\n", "'''convert geospatial coordinates to pixels'''\n", "def world2screen(bbox, w, h, x, y):\n", " minx,miny,maxx,maxy = bbox\n", " xdist = maxx - minx\n", " ydist = maxy - miny\n", " xratio = w/xdist\n", " yratio = h/ydist\n", " px = int(w - ((maxx - x) * xratio))\n", " py = int((maxy - y) * yratio)\n", " return (px,py)\n", "\n", "# Open our shapefile\n", "inShp = shapefile.Reader(\"GIS_CensusTract_poly\")\n", "iwidth = 600\n", "iheight = 400\n", "\n", "# PIL Image\n", "img = Image.new(\"RGB\", (iwidth,iheight), (255,255,255)) \n", "\n", "# PIL Draw module for polygon fills\n", "draw = ImageDraw.Draw(img)\n", "\n", "# Get the population AND area index\n", "pop_index = None\n", "area_index = None\n", "\n", "# Shade the census tracts\n", "for i,f in enumerate(inShp.fields):\n", " if f[0] == \"POPULAT11\":\n", " # Account for deletion flag\n", " pop_index = i-1\n", " elif f[0] == \"AREASQKM\":\n", " area_index = i-1\n", "\n", "# Draw the polygons\n", "for sr in inShp.shapeRecords():\n", " density = sr.record[pop_index]/sr.record[area_index]\n", " weight = min(math.sqrt(density/80.0), 1.0) * 50 \n", " R = int(205 - weight)\n", " G = int(215 - weight)\n", " B = int(245 - weight)\n", " pixels = []\n", " for x,y in sr.shape.points:\n", " (px,py) = world2screen(inShp.bbox, iwidth, iheight, x, y)\n", " pixels.append((px,py))\n", " draw.polygon(pixels, outline=(255,255,255), fill=(R,G,B))\n", " \n", "# Save the image to file\n", "img.save(\"choropleth.png\")" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Shapefile contains 1 layers\n", "\n", "Layer 0 has spatial reference +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs \n", "Layer 0 has 955 features\n", "\n", "Feature 0 has name Hilo, HI\n", "Feature 1 has name Bozeman, MT\n", "Feature 2 has name Magnolia, AR\n", "Feature 3 has name Kingston, NY\n", "Feature 4 has name Astoria, OR\n", "Feature 5 has name Kennewick-Pasco-Richland, WA\n", "Feature 6 has name Fayetteville-Springdale-Rogers, AR-MO\n", "Feature 7 has name Seattle-Tacoma-Bellevue, WA\n", "Feature 8 has name Ellensburg, WA\n", "Feature 9 has name Pierre, SD\n" ] } ], "source": [ "from osgeo import ogr\n", "\n", "shapefile_name = \"tl_2012_us_cbsa.shp\"\n", "shapefile = ogr.Open(shapefile_name)\n", "numLayers = shapefile.GetLayerCount()\n", "\n", "print (\"Shapefile contains %d layers\\n\" % numLayers)\n", "\n", "for layerNum in range(numLayers):\n", " layer = shapefile.GetLayer(layerNum)\n", " spatialRef = layer.GetSpatialRef().ExportToProj4()\n", " numFeatures = layer.GetFeatureCount()\n", " print (\"Layer %d has spatial reference %s\" % (layerNum, spatialRef))\n", " print (\"Layer %d has %d features\\n\" % (layerNum, numFeatures))\n", "\n", "for featureNum in range(10):\n", " feature = layer.GetFeature(featureNum) \n", " featureName = feature.GetField(\"NAME\")\n", " print (\"Feature %d has name %s\" % (featureNum, featureName))\n", "\n" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Feature 2 has the following attributes:\n", " CSAFP = None\n", " CBSAFP = 31620\n", " GEOID = 31620\n", " NAME = Magnolia, AR\n", " NAMELSAD = Magnolia, AR Micro Area\n", " LSAD = M2\n", " MEMI = 2\n", " MTFCC = G3110\n", " ALAND = 1984070931\n", " AWATER = 1809509\n", " INTPTLAT = +33.2230377\n", " INTPTLON = -093.2328433\n", "\n", "Feature's geometry data consists of a POLYGON\n" ] } ], "source": [ "from osgeo import ogr\n", " \n", "shapefile_name = \"tl_2012_us_cbsa.shp\"\n", "shapefile = ogr.Open(shapefile_name)\n", "layer = shapefile.GetLayer(0)\n", "feature = layer.GetFeature(2)\n", " \n", "print (\"Feature 2 has the following attributes:\")\n", " \n", "attributes = feature.items()\n", " \n", "for key,value in attributes.items():\n", " print (\" %s = %s\" % (key, value))\n", " \n", "geometry = feature.GetGeometryRef()\n", "geometryName = geometry.GetGeometryName()\n", "\n", "print ()\n", "print (\"Feature's geometry data consists of a %s\" % geometryName)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "s = []\n", "s.append(\" \" * 0)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON containing:\n" ] } ], "source": [ "s = []\n", "s.append(\" \" * 0)\n", "s.append(geometry.GetGeometryName())\n", "if geometry.GetPointCount() > 0:\n", " s.append(\" with %d data points\" % geometry.GetPointCount())\n", "if geometry.GetGeometryCount() > 0:\n", " s.append(\" containing:\")\n", " \n", "print (\"\".join(s))" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON containing:\n", " LINEARRING with 9193 data points\n" ] } ], "source": [ "from osgeo import ogr\n", "\n", "def analyzeGeometry(geometry, indent=0):\n", " s = []\n", " s.append(\" \" * indent)\n", " s.append(geometry.GetGeometryName())\n", " if geometry.GetPointCount() > 0:\n", " s.append(\" with %d data points\" % geometry.GetPointCount())\n", " if geometry.GetGeometryCount() > 0:\n", " s.append(\" containing:\")\n", "\n", " print (\"\".join(s))\n", "\n", " for i in range(geometry.GetGeometryCount()):\n", " analyzeGeometry(geometry.GetGeometryRef(i), indent+1)\n", "\n", "shapefile_name = \"tl_2012_us_cbsa.shp\"\n", "shapefile = osgeo.ogr.Open(shapefile_name)\n", "layer = shapefile.GetLayer(0)\n", "feature = layer.GetFeature(55)\n", "geometry = feature.GetGeometryRef()\n", "\n", "analyzeGeometry(geometry)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 1.0 First\n", "3.0 1.0 Second\n", "4.0 3.0 Third\n", "2.0 2.0 Fourth\n", "0.0 0.0 Appended\n" ] } ], "source": [ "# Examine a shapefile with pyshp \n", "import shapefile\n", " \n", "shp = shapefile.Reader(\"point\")\n", "for feature in shp.shapeRecords():\n", " point = feature.shape.points[0]\n", " rec = feature.record[0]\n", " print (point[0], point[1], rec)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Northernmost point is (-93.9197, 33.0195)\n", "Southernmost point is (-93.8752, 31.8442)\n" ] } ], "source": [ "from osgeo import ogr\n", "\n", "def findPoints(geometry, results):\n", " for i in range(geometry.GetPointCount()):\n", " x,y,z = geometry.GetPoint(i)\n", " if results['north'] == None or results['north'][1] < y:\n", " results['north'] = (x,y)\n", " if results['south'] == None or results['south'][1] > y:\n", " results['south'] = (x,y)\n", "\n", " for i in range(geometry.GetGeometryCount()):\n", " findPoints(geometry.GetGeometryRef(i), results)\n", "\n", "shapefile_name = \"tl_2012_us_cbsa.shp\"\n", "shapefile = ogr.Open(shapefile_name)\n", "layer = shapefile.GetLayer(0)\n", "feature = layer.GetFeature(55)\n", "geometry = feature.GetGeometryRef()\n", "\n", "results = {'north' : None,\n", " 'south' : None}\n", "\n", "findPoints(geometry, results)\n", "\n", "print (\"Northernmost point is (%0.4f, %0.4f)\" % results['north'])\n", "print (\"Southernmost point is (%0.4f, %0.4f)\" % results['south'])\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Sources : this material was adopted from the following sources\n", "- https://github.com/ycrc/Data-Analysis-with-Python\n", "- https://github.com/fpl/geotutorial_basic\n", "- Python GeoSpatial Development: https://s3.amazonaws.com/academia.edu.documents/32168914/Python_Geospatial_Development_%282010%29.pdf?AWSAccessKeyId=AKIAIWOWYYGZ2Y53UL3A&Expires=1534884658&Signature=3bXh%2BzqzQJd%2BBaLgYNVxRmq2K8c%3D&response-content-disposition=inline%3B%20filename%3DPython_Geospatial_Development_2010.pdf" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Slideshow", "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": 2 }