Created
December 14, 2016 21:33
-
-
Save mbjoseph/9934454330b5715db9529502cf415f72 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## pyDEM failure on a SRTM 30m DEM \n", | |
"\n", | |
"This is a minimal example to demonstrate a use case where slope and aspect calculations seem to fail. " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"import elevation\n", | |
"import matplotlib.pyplot as plt\n", | |
"import matplotlib\n", | |
"import numpy as np\n", | |
"from osgeo import gdal\n", | |
"from pydem.dem_processing import DEMProcessor\n", | |
"\n", | |
"%matplotlib inline\n", | |
"matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Using the `elevation` command line interface to get a DEM: " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"make: Nothing to be done for `download'.\n", | |
"make: Nothing to be done for `all'.\n", | |
"gdal_translate -q -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 -co PREDICTOR=2 -projwin -122.4 41.5 -122.1 41.2 SRTM1.vrt /Users/majo3748/Documents/earthlab/tutorials/Shasta-30m-DEM.tif\n" | |
] | |
} | |
], | |
"source": [ | |
"!eio clip -o Shasta-30m-DEM.tif --bounds -122.4 41.2 -122.1 41.5 " | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Print projection information and verify WGS 84." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]\n" | |
] | |
} | |
], | |
"source": [ | |
"filename = \"Shasta-30m-DEM.tif\"\n", | |
"gdal_data = gdal.Open(filename)\n", | |
"prj = gdal_data.GetProjection()\n", | |
"print(prj)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Now, computing slope and aspect with `pyDEM`:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"starting slope/direction calculation for chunk 1 [0:516, 0:516]\n" | |
] | |
}, | |
{ | |
"name": "stderr", | |
"output_type": "stream", | |
"text": [ | |
"/Users/majo3748/anaconda/envs/py27/lib/python2.7/site-packages/pydem/dem_processing.py:563: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.\n", | |
"Check the NumPy 1.11 release notes for more information.\n", | |
" self.data[data.mask] = FILL_VALUE\n" | |
] | |
}, | |
{ | |
"ename": "TypeError", | |
"evalue": "Cannot cast ufunc add output from dtype('float64') to dtype('bool') with casting rule 'same_kind'", | |
"output_type": "error", | |
"traceback": [ | |
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", | |
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", | |
"\u001b[0;32m<ipython-input-4-9def4a7301a2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mprocessor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDEMProcessor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mmag\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maspect\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprocessor\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalc_slopes_directions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", | |
"\u001b[0;32m/Users/majo3748/anaconda/envs/py27/lib/python2.7/site-packages/pydem/dem_processing.pyc\u001b[0m in \u001b[0;36mcalc_slopes_directions\u001b[0;34m(self, plotflag)\u001b[0m\n\u001b[1;32m 856\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mte\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mbe\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 857\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdY\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mte\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mbe\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmag\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 858\u001b[0;31m direction)\n\u001b[0m\u001b[1;32m 859\u001b[0m \u001b[0mdirection\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mflats\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFLAT_ID\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 860\u001b[0m \u001b[0mmag\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mflats\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFLAT_ID\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", | |
"\u001b[0;32m/Users/majo3748/anaconda/envs/py27/lib/python2.7/site-packages/pydem/dem_processing.pyc\u001b[0m in \u001b[0;36m_find_flats_edges\u001b[0;34m(self, data, dX, dY, mag, direction)\u001b[0m\n\u001b[1;32m 939\u001b[0m \u001b[0mids_tmp2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0melev_flat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 940\u001b[0m \u001b[0mflats\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mids_tmp2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 941\u001b[0;31m \u001b[0;34m+=\u001b[0m \u001b[0mFLATS_KERNEL3\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miii\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjjj\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 942\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 943\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mflats\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", | |
"\u001b[0;31mTypeError\u001b[0m: Cannot cast ufunc add output from dtype('float64') to dtype('bool') with casting rule 'same_kind'" | |
] | |
} | |
], | |
"source": [ | |
"processor = DEMProcessor(filename)\n", | |
"mag, aspect = processor.calc_slopes_directions()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"plt.imshow(np.log(mag), cmap=\"magma\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"plt.imshow(aspect, cmap = \"viridis\")" | |
] | |
} | |
], | |
"metadata": { | |
"anaconda-cloud": {}, | |
"kernelspec": { | |
"display_name": "Python [py27]", | |
"language": "python", | |
"name": "Python [py27]" | |
}, | |
"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.11" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 0 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment