tag:blogger.com,1999:blog-20445554334063981562024-03-17T22:03:45.882-05:00GeospatialPython.comPoints, Polylines, Polygons, Pixels, Python!Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.comBlogger71125tag:blogger.com,1999:blog-2044555433406398156.post-50898207983301952252023-10-30T14:09:00.004-05:002023-10-30T14:20:49.410-05:00Thank you for making "Learning Geospatial Analysis with Python" #1 on Amazon!<p> <span color="rgba(0, 0, 0, 0.9)" face="-apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif" style="background-color: white; font-size: 14px;">🚀 Exciting News!! 📚 The 10-year anniversary edition of 🌎 "Learning Geospatial Analysis with Python - 4th Ed," 🐍 is now available on Amazon and it's starting off on fire 🔥 as #1 in the GIS New Release category! 🥇</span></p><span color="rgba(0, 0, 0, 0.9)" face="-apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif" style="background-color: white; font-size: 14px;">Thank you to all of my fellow geography-obsessed and Python-loving friends out there for making this book possible. In this edition, I'll help you up your geospatial game with package management in Anaconda and new data types including bathymetric point clouds so you can explore the ocean. Even more exciting I'll show you how to speed up your workflow by using ChatGPT as a tireless geospatial programming assistant.</span><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><span color="rgba(0, 0, 0, 0.9)" face="-apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif" style="background-color: white; font-size: 14px;">I bring you up to speed on the entire history of geospatial analysis and then lead you through concepts and major software tools, and then teach you simple geospatial algorithms that will serve as building blocks for the complete programs that come later in the book.</span><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><span color="rgba(0, 0, 0, 0.9)" face="-apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif" style="background-color: white; font-size: 14px;">These examples are designed to teach people new to Python and geospatial analysis, but are built to empower experienced practitioners as well with code that you will use everyday for common GIS challenges.</span><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><span color="rgba(0, 0, 0, 0.9)" face="-apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif" style="background-color: white; font-size: 14px;">Grab 🖐 your copy 📖 today and start changing the world by analyzing it using one of the most popular programming languages ever created!</span><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><br style="background-color: white; box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit;" /><a data-attribute-index="0" href="https://a.co/d/atLzCEc" style="background-color: white; border: var(--artdeco-reset-link-border-zero); box-sizing: inherit; font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; font-weight: var(--font-weight-bold); line-height: inherit; margin: var(--artdeco-reset-base-margin-zero); overflow-wrap: normal; padding: var(--artdeco-reset-base-padding-zero); position: relative; text-decoration: var(--artdeco-reset-link-text-decoration-none); touch-action: manipulation; vertical-align: var(--artdeco-reset-base-vertical-align-baseline); word-break: normal;" target="_self">https://a.co/d/atLzCEc</a><div><br /></div><div class="separator" style="clear: both; text-align: left;"><a href="https://a.co/d/atLzCEc" style="margin-left: 1em; margin-right: 1em;" target="_blank"><img border="0" data-original-height="749" data-original-width="1774" height="235" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhDVhlNoYbby39xeBnS5Z-y-LlBPJFXsaQAuZUiOCqzqhJUS3LbWG83b52ojDrtKReeC47QeBctQxUuF-kOaNxGfVLIfOEM0esj61xqtR_DY8le8TlkbQwJOeUdsXzdHZDeUnR_uoG3oHNKKPBVHMx5TSlZ3bfhlC21nSL4tKXJUzMSaNhqHXP7FpHySDw/w558-h235/IMG_5173.jpg" width="558" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://a.co/d/atLzCEc" style="margin-left: 1em; margin-right: 1em;" target="_blank"><img border="0" data-original-height="658" data-original-width="1699" height="201" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTV8Q6MMEV6fyyLwYphsc1sOTO_65aIgtY4nUyU75RzHXmDJ6nzMVDnS7izWbWom44W68Q0bucBXygX23RzOqX2hvZc7L1TGk9v9S6OY6YcCSdd1PZgBrAQsYIGjQ7BEavCiL2nLhS-6dcWyLcZn6DkPCuMlje3BVmf8gDy6AeXwXEjLGnh6KS3KDNZ_4/w517-h201/IMG_5174.jpg" width="517" /></a></div><br /><div><br /></div>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-51191604110095388402023-10-26T20:50:00.003-05:002023-10-26T20:50:28.999-05:00Learning Geospatial Analysis with Python: 10-year Anniversary Edition!<p><span style="background-color: white; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px;">I'm happy to announce the 10-year anniversary of the 1st edition of my book, "Learning Geospatial Analysis with Python". There's a lot of new information in this 4th edition. I've always tried to provide pure Python solutions, but this edition introduces the amazing Anaconda platform which makes installing complex geospatial libraries a breeze. I've also introduced my first example of dealing with hydrographic data as the "Blue Economy" continues to rise. And most importantly, I provide examples of using ChatGPT to assist in quickly creating useful geospatial analysis Python scripts as we move into the age of AI. Available on Amazon:</span><span class="white-space-pre" style="background-color: white; border: var(--artdeco-reset-base-border-zero); box-sizing: inherit; color: rgba(0, 0, 0, 0.9); font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; line-height: inherit !important; margin: var(--artdeco-reset-base-margin-zero); outline: var(--artdeco-reset-base-outline-zero); padding: var(--artdeco-reset-base-padding-zero); vertical-align: var(--artdeco-reset-base-vertical-align-baseline); white-space: pre !important;"> </span><a class="app-aware-link " data-test-app-aware-link="" href="https://lnkd.in/gBuvbvXD" style="background-color: white; border: var(--artdeco-reset-link-border-zero); box-sizing: inherit; font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; font-weight: var(--font-weight-bold); line-height: inherit !important; margin: var(--artdeco-reset-base-margin-zero); overflow-wrap: normal; padding: var(--artdeco-reset-base-padding-zero); position: relative; text-decoration: var(--artdeco-reset-link-text-decoration-none); touch-action: manipulation; vertical-align: var(--artdeco-reset-base-vertical-align-baseline); word-break: normal;" target="_self">https://lnkd.in/gBuvbvXD</a></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg8rzOYz-AWkwcP9fiXMUayFZakQ19qcg58gz52KpKi_S4Uv8t1u5weZQ8LLDR5voWT-_LZsKCJk3tbCoQzGCM7sZEhjCYpvlCXtIUqJ-A6DYtwlQoJ2lEzDOFxHFwAqh0GxihSowTDIcZuJ1-wA9rgJeJX4pid08njA0yZn16bqq9WOA_TY-aI2lglHZw/s1500/LGAWP4thEdCover.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1500" data-original-width="1216" height="563" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg8rzOYz-AWkwcP9fiXMUayFZakQ19qcg58gz52KpKi_S4Uv8t1u5weZQ8LLDR5voWT-_LZsKCJk3tbCoQzGCM7sZEhjCYpvlCXtIUqJ-A6DYtwlQoJ2lEzDOFxHFwAqh0GxihSowTDIcZuJ1-wA9rgJeJX4pid08njA0yZn16bqq9WOA_TY-aI2lglHZw/w456-h563/LGAWP4thEdCover.jpg" width="456" /></a></div><br /><p><br /></p>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-49031084076809607382021-08-20T06:47:00.003-05:002021-10-18T09:51:32.835-05:00Connecting MultiPolygon Edges to Lines<div><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiupZlBagJekGG80qCYey0BxO0divWjqGTBj3FYHP7ODVulJmRvY0CRrntD6ERCUBIWF8E7cfgpkySKV8qvAyOYZDL6Gmbk4sHN94m9KbKMeUwfk4DswqDtsxyrE6hhLnFsJxhyovt0-NI/s1031/1629223501030.gif" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" data-original-height="878" data-original-width="1031" height="274" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiupZlBagJekGG80qCYey0BxO0divWjqGTBj3FYHP7ODVulJmRvY0CRrntD6ERCUBIWF8E7cfgpkySKV8qvAyOYZDL6Gmbk4sHN94m9KbKMeUwfk4DswqDtsxyrE6hhLnFsJxhyovt0-NI/w320-h274/1629223501030.gif" width="320" /></a>This post is another <a href="https://www.linkedin.com/posts/spatial-thoughts_spatialanalysis-activity-6833458650091274240-yOI3" target="_blank">Spatial Thoughts Academy Weekly Challenge</a> solutions. The challenge is to find the edge of the polygon in a set of building footprints whose center point is closest to a street and then connect it with the nearest point on the street. Shapely and Fiona were of course up to the challenge. The basic algorithm is to grab the exterior-most ring of the building footprints, break out the edges, determine the midpoint of each edge, measure the distance to each road until the shortest distance is found, then grab the nearest point on the road, and finally draw a line between the edge midpoint and the nearest road point. Shapely's nearest_points() method is really versatile here because you can throw any two geometries at it even if they are different types. Fiona is used to read and write geopackages which I'm really liking lately as more an more people are using them finally. Here's the code:</div><div><br /></div><script src="https://gist.github.com/GeospatialPython/d65b98ed1122621e185733f4780eead9.js"></script>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-52260401579059428392021-08-12T05:11:00.002-05:002021-08-12T05:12:40.957-05:00Extracting MultiPolygons with Holes using Shapely and Fiona<div class="separator"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgvTOQb8651MuSogjO2BlPrYOzX96-dNJywMzdafQBWezDXK0loquNf9Y3hMZl9WfzV3_3wdYm6Kr6oGeB4JXfcTqI8vcGCeCuOwp4X-iUtLc3nxydFAVP7lVC0jNecZ2NJ7SMIU0uzz_0/s800/BuildingFootprints.gif" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em; text-align: center;"><img border="0" data-original-height="558" data-original-width="800" height="279" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgvTOQb8651MuSogjO2BlPrYOzX96-dNJywMzdafQBWezDXK0loquNf9Y3hMZl9WfzV3_3wdYm6Kr6oGeB4JXfcTqI8vcGCeCuOwp4X-iUtLc3nxydFAVP7lVC0jNecZ2NJ7SMIU0uzz_0/w400-h279/BuildingFootprints.gif" width="400" /></a></div><a href="https://spatialthoughts.com">SpatialThoughts.com</a> recently <a href="https://www.linkedin.com/posts/spatial-thoughts_spatialanalysis-activity-6830545160842158081--FMH" target="_blank">posted a challenge</a> on LinkedIn to extract only building footprints with holes from a city-wide dataset. <div><br /></div><div>What made the challenge interesting is the building footprints are MultiPolygons instead of simple Polygons which stumped a few people's attempts. But this challenge is very straight-forward Shapely to select the polygons and Fiona to read and write the data.</div><div><br /></div><div>The correct feature count if you extract all of the building footprints with holes or polygon inner rings is 45. The challenge author later demonstrated how to do it in QGIS using a series of expressions. I think I had the only Python solution. One guy posted a successful answer in R.</div><div><br /></div><div>The challenge data is a city extracted from OSM: <a href="https://lnkd.in/g7dVqzMx">https://lnkd.in/g7dVqzMx</a></div><div><br /></div><div>Here is my solution:</div><div></div><div><br /></div>
<script src="https://gist.github.com/GeospatialPython/a7500c04741a568bcd0a6ee3e1c6902c.js"></script>Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-27953529655536950882019-10-10T05:38:00.002-05:002019-10-10T05:38:50.473-05:00Learning Geospatial Analysis with Python, 3rd Ed.<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhS-u_7dDftyZ0oGWeuSJneuhGDLcIfHX2Shxkn_tY6_pMjMxUNKHLBv_imtdmkjJCdC76fyJ9viGa1diQMxGaQnUwOns_C5muMB7K9hMIgK9Qnhy1mHYT6QnnmiGOnwZ05Uh_wx65y2gU/s1600/9781789959277-original.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" data-original-height="1000" data-original-width="810" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhS-u_7dDftyZ0oGWeuSJneuhGDLcIfHX2Shxkn_tY6_pMjMxUNKHLBv_imtdmkjJCdC76fyJ9viGa1diQMxGaQnUwOns_C5muMB7K9hMIgK9Qnhy1mHYT6QnnmiGOnwZ05Uh_wx65y2gU/s320/9781789959277-original.jpg" width="259" /></a></div>
<span class="ember-view" id="ember3495" style="border: 0px; box-sizing: inherit; color: rgba(0 , 0 , 0 , 0.901961); font-family: , , "blinkmacsystemfont" , "segoe ui" , "roboto" , "helvetica neue" , "fira sans" , "ubuntu" , "oxygen" , "oxygen sans" , "cantarell" , "droid sans" , "apple color emoji" , "segoe ui emoji" , "segoe ui symbol" , "lucida grande" , "helvetica" , "arial" , sans-serif; font-size: 14px; line-height: inherit; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline; white-space: pre-wrap;">Third Edition is on the shelves! Geospatial concepts, Geo-python universe, and pound-for-pound still the most pure-python and minimal-dependency examples you’ll find anywhere so somebody somewhere out there will still be able to do the math. I'm giving away a few eBook versions if you're willing to write a review on Amazon. DM me @SpatialPython with your name and email! </span><a class="feed-shared-text-view__hyperlink ember-view" href="https://bit.ly/35gm0xl" id="ember3499" rel="noopener noreferrer" style="border: 0px; box-sizing: inherit; color: #0073b1; font-family: -apple-system, system-ui, BlinkMacSystemFont, "Segoe UI", Roboto, "Helvetica Neue", "Fira Sans", Ubuntu, Oxygen, "Oxygen Sans", Cantarell, "Droid Sans", "Apple Color Emoji", "Segoe UI Emoji", "Segoe UI Symbol", "Lucida Grande", Helvetica, Arial, sans-serif; font-size: 14px; font-weight: 600; line-height: inherit !important; margin: 0px; padding: 0px; position: relative; text-decoration: none; touch-action: manipulation; vertical-align: baseline; white-space: pre-wrap;" target="_blank">https://bit.ly/35gm0xl</a>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-59632076941039700582019-01-12T10:34:00.000-06:002019-01-12T10:34:03.433-06:00PotreeConverter Mac Binary<a href="http://www.potree.org/" target="_blank">Potree</a> is the amazing javascript WebGL library that can effortlessly display multi-million-point lidar <div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgeF2so4nwd8Q3xVDjNRdHyh-u3l-DGcrs6IwY2DBZw8XevEo07pEUuWu876QWXroK6YKtuNKbTlZdr8MxPIcguNcsqiNWwOM8Xm90TyQUhTQ9_4n4NbdwTZda0l-899QgRuzVNMuhIsTo/s1600/Screen+Shot+2019-01-12+at+10.05.14+AM.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" data-original-height="1000" data-original-width="1600" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgeF2so4nwd8Q3xVDjNRdHyh-u3l-DGcrs6IwY2DBZw8XevEo07pEUuWu876QWXroK6YKtuNKbTlZdr8MxPIcguNcsqiNWwOM8Xm90TyQUhTQ9_4n4NbdwTZda0l-899QgRuzVNMuhIsTo/s320/Screen+Shot+2019-01-12+at+10.05.14+AM.png" width="320" /></a></div>
point clouds in a browser using a static web page. In order to do that, you have to use <a href="https://github.com/potree/PotreeConverter" target="_blank">PotreeConverter</a> which creates an efficient octree of your LAS file. But the MacOS built-in compilers with XCode are usually several generations behind and don't like the PotreeConverter source. I was finally able to get it to compile so I created a <a href="https://github.com/GeospatialPython/MacOS-PotreeConverter/" target="_blank">GitHub repository</a> with the binary and the compile steps. You still need to download the "resources" directory from the original PotreeConverter repository which contains the viewer webpage template. You may also have to deal with issue <a href="https://github.com/potree/PotreeConverter/issues/281">https://github.com/potree/PotreeConverter/issues/281</a> but it's an easy fix. I have no idea how portable the binary is, but hopefully the compile steps will save time for others. <div class="separator" style="clear: both; text-align: center;">
</div>
Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-43197336865739833752017-11-16T16:20:00.001-06:002017-11-16T16:21:06.345-06:00Creating MultiPoint Shapefiles with PyShpPyshp let's you create any type of shapefile. Normally a point shapefile has one point per record. But<br />
<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: right; margin-left: 1em; text-align: right;"><tbody>
<tr><td style="text-align: center;"><a href="http://www.vancouvermom.ca/wp-content/uploads/2015/04/ball-pit_Lonsdale2.jpg" imageanchor="1" style="clear: right; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" src="http://www.vancouvermom.ca/wp-content/uploads/2015/04/ball-pit_Lonsdale2.jpg" data-original-height="535" data-original-width="800" height="132" width="200" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Photo credit: <a href="http://vancouvermom.ca/">VancouverMom.ca</a></td></tr>
</tbody></table>
you can also have a group of points tied to a single record in a Multipoint shapefile. To create a Multipoint shapefile, you just use the "poly" method in the Writer. The poly method isn't just for polygons. It can also create polylines and polypoints (i.e. MultiPoints)!<br />
<br />
Here's a simple Multipoint shapefile:<br />
<br />
<br />
<pre class="prettyprint">import shapefile
# Create our writer as a multipoint shapefile
w = shapefile.Writer(shapefile.MULTIPOINT)
# We'll create a single dbf field
w.field("NAME","C",40)
# Create a single-part, multi-point shape
# by declaring the shape
# type after the parts/points list
w.poly([[[3,4],[5,6],[7,8],[9,8]]], shapeType=shapefile.MULTIPOINT)
# Create a record for this feature
w.record("Group1")
# Save the multipoint shapefile
w.save("mpoint")</pre>
Repeat the poly and record steps to add additional shapes. Add another nested list of points in the first poly method argument to add more parts to the same recordJoel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-88666900327999347882017-03-14T15:12:00.002-05:002017-03-14T15:12:24.345-05:00QGIS Python Programming Cookbook - SECOND EDITION<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4DNPdVl7UQnn_Cqbehp1E68qcYwNJdRmS5e1iq6NQq0z3Fsvqas0RsYQmEBA-F4KGP8q5EHub9X6t-nVURnLp3fOgm74_breWCOhxbNXhg3bwnAOJ6LKniWKEqi53u0-_ng2u7KbfBR4/s1600/book.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4DNPdVl7UQnn_Cqbehp1E68qcYwNJdRmS5e1iq6NQq0z3Fsvqas0RsYQmEBA-F4KGP8q5EHub9X6t-nVURnLp3fOgm74_breWCOhxbNXhg3bwnAOJ6LKniWKEqi53u0-_ng2u7KbfBR4/s320/book.png" width="259" /></a>The second edition of the "<a href="https://www.packtpub.com/application-development/qgis-python-programming-cookbook-second-edition" target="_blank">QGIS Python Programming Cookbook</a>" is out today from Packt Publishing! And after my second publishing experience writing about QGIS, I can enthusiastically say it is one of the greatest open source projects, GIS software, and Python APIs out there. It really sets a standard for quality software development. GIS is a tool that helps people make better decisions to make the world a better place, and QGIS puts that potential in more people's hands than any other GIS software.<br />
<br />
As with all of my books, I've tried to walk a very tricky line. I created a QGIS Python book reference book that will stand on its own if it's the only book you have. But I also try to create content that is very different from the official documentation and book, as well as complement other third-party books out there. QGIS is rich enough, <br />
and updated and released with such frequency, that I don't think one book can do it justice. While many operations in QGIS are straight-forward and easy to understand, some things can be quite complex. For example, the Map Composer, raster operations, tinkering with various types of settings, and working with plug-ins through Python can be more difficult to understand. And having more than one example from different sources can be invaluable in software development. So I carefully picked examples and explanations that I felt were not easily found on the Internet for the most part, lacked documentation, or were complex enough that the over 170 examples should enrich your knowledge of QGIS regardless of you're starting from. At the very least, I hope it encourages more people to use my favorite GIS package. Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-31055786517267216942016-10-16T21:56:00.001-05:002016-10-16T22:03:11.304-05:00Cutting a Line Shapefile with a Point Shapefile using Shapely<a href="https://github.com/Toblerity/Shapely" target="_blank">Shapely</a> works with geometric classes. Sometimes you want those geometries to come from<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbPLZrTWOlopHe6sw6ckJiVTBOJZxGWEc-ZqFbGXBqH0f-PMtdv7oiatzOZgXAkITyDDm3-aNlya0aA8df1ACE6TjBGNzCRhZgsKn0CV55LizrHKqUZS8ogDdV-NLh2-fPFvUGWpk3Ls4/s1600/Untitled.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="177" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbPLZrTWOlopHe6sw6ckJiVTBOJZxGWEc-ZqFbGXBqH0f-PMtdv7oiatzOZgXAkITyDDm3-aNlya0aA8df1ACE6TjBGNzCRhZgsKn0CV55LizrHKqUZS8ogDdV-NLh2-fPFvUGWpk3Ls4/s200/Untitled.png" width="200" /></a></div>
shapefiles. That's where PyShp comes in. In this example on GIS Stack Exchange, we take a line shapefile and a point shapefile with points on or near the line segments, and split those lines at the location of the nearest point. At a minimum, the point must be inside the envelope of the line segment:<br />
<br />
<a href="http://gis.stackexchange.com/questions/214344/split-polyline-shapefile-using-python-with-point-shapefile/214432#214432">http://gis.stackexchange.com/questions/214344/split-polyline-shapefile-using-python-with-point-shapefile/214432#214432</a><br />
<br />Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-6319277031879949122016-10-12T14:40:00.001-05:002016-10-12T14:40:40.646-05:00Reproject a Polygon Shapefile using PyShp and PyProj<span style="font-family: Arial, Helvetica, sans-serif;">Via the blog "Geospatiality": <span style="background-color: white; color: #373737; font-size: 15px;">In this post I will use the </span><em style="background-color: white; border: 0px; color: #373737; font-size: 15px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;"><strong style="border: 0px; font-style: inherit; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">PyShp</strong> </em><span style="background-color: white; color: #373737; font-size: 15px;">library along with the </span><em style="background-color: white; border: 0px; color: #373737; font-size: 15px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;"><strong style="border: 0px; font-style: inherit; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">PyProj</strong> </em><span style="background-color: white; color: #373737; font-size: 15px;">library to </span></span><br />
<div class="separator" style="clear: both; text-align: center;">
<span style="font-family: Arial, Helvetica, sans-serif;"><a href="https://geobitz.files.wordpress.com/2016/01/itm_to_wgs84.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="135" src="https://geobitz.files.wordpress.com/2016/01/itm_to_wgs84.png" width="320" /></a></span></div>
<span style="font-family: Arial, Helvetica, sans-serif;">reproject the <em style="background-color: white; border: 0px; color: #373737; font-size: 15px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">local authority boundarie</em><span style="background-color: white; color: #373737; font-size: 15px;">s of Ireland, in Shapefile format, from </span><strong style="background-color: white; border: 0px; color: #373737; font-size: 15px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">Irish Transverse Mercator to WGS 84 </strong><span style="background-color: white; color: #373737; font-size: 15px;">using </span><em style="background-color: white; border: 0px; color: #373737; font-size: 15px; margin: 0px; outline: 0px; padding: 0px; vertical-align: baseline;">Python</em><span style="background-color: white; color: #373737; font-size: 15px;">.</span></span><br />
<span style="background-color: white; color: #373737; font-size: 15px;"><span style="font-family: Arial, Helvetica, sans-serif;"><br /></span></span>
<span style="background-color: white; font-size: 15px;"><span style="color: #373737; font-family: Arial, Helvetica, sans-serif;"><a href="https://glenbambrick.com/2016/01/24/reproject-shapefile/">https://glenbambrick.com/2016/01/24/reproject-shapefile/</a></span></span>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-36206832047001009962016-10-04T16:55:00.002-05:002016-10-04T16:58:26.616-05:003D Multipatch Shapefile in Pure PythonVia GIS StackExchange...<br />
<br />
<a href="http://gis.stackexchange.com/questions/212821/pyshp-create-multipatch-cube/">http://gis.stackexchange.com/questions/212821/pyshp-create-multipatch-cube/</a><br />
<br />
<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: left; margin-right: 1em; text-align: left;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSxbQCHX1Zlf0aPovf_Cpv1GMy6APmhq5zIjqcKezopyIbI-fJsJrtpU7TAtnVQcuoqxA2AT700_F8nNd_1M8x-5r4hKSgwrf1caPJwlJLB6Sw1cZ_63V-dfB5JpsRL0PH4Jf6q3rnNwM/s1600/cube.png" imageanchor="1" style="clear: left; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="237" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSxbQCHX1Zlf0aPovf_Cpv1GMy6APmhq5zIjqcKezopyIbI-fJsJrtpU7TAtnVQcuoqxA2AT700_F8nNd_1M8x-5r4hKSgwrf1caPJwlJLB6Sw1cZ_63V-dfB5JpsRL0PH4Jf6q3rnNwM/s400/cube.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><b>3D cube shown in Esri ArcScene generated using PyShp <br />and the obscure "MultiPatch" shapefile type in pure Python</b></td></tr>
</tbody></table>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSxbQCHX1Zlf0aPovf_Cpv1GMy6APmhq5zIjqcKezopyIbI-fJsJrtpU7TAtnVQcuoqxA2AT700_F8nNd_1M8x-5r4hKSgwrf1caPJwlJLB6Sw1cZ_63V-dfB5JpsRL0PH4Jf6q3rnNwM/s1600/cube.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"></a><br /></div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSxbQCHX1Zlf0aPovf_Cpv1GMy6APmhq5zIjqcKezopyIbI-fJsJrtpU7TAtnVQcuoqxA2AT700_F8nNd_1M8x-5r4hKSgwrf1caPJwlJLB6Sw1cZ_63V-dfB5JpsRL0PH4Jf6q3rnNwM/s1600/cube.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"></a>Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-66306991883659930352016-08-07T18:24:00.001-05:002016-08-08T12:36:19.535-05:00Pure Python: Loading Shapefiles into PostGIS<div class="separator" style="clear: both; text-align: center;">
</div>
<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: right; margin-left: 1em; text-align: right;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgTE3E2g3QCFcfleQWYpy0wsjor6BvAw623ToOLRplVsUU7Amz_ICis6PI8fEtOCP2_a3Q1TkniATWk9OmvdwKz9C2zITW04U4-EpdnzMuV33QN3md0oKx252YBhjxCl_GTadAOMEdrNOo/s1600/postgis.png" imageanchor="1" style="clear: right; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="305" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgTE3E2g3QCFcfleQWYpy0wsjor6BvAw623ToOLRplVsUU7Amz_ICis6PI8fEtOCP2_a3Q1TkniATWk9OmvdwKz9C2zITW04U4-EpdnzMuV33QN3md0oKx252YBhjxCl_GTadAOMEdrNOo/s320/postgis.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><div style="text-align: left;">
<b>In addition to being the best spatial database <br />engine </b><b style="font-size: 12.8px;"><b>available, PostGIS can be manipulated <br />using pure Python </b></b><b style="font-size: 12.8px;"><b>without any C-library <br />dependencies.</b></b></div>
<b></b><br />
<b>
</b></td></tr>
</tbody></table>
<a href="http://postgis.net/" target="_blank">PostGIS</a> is the spatial database engine for the open source <a href="https://www.postgresql.org/" target="_blank">PostgreSQL</a> database. For a decade we<a href="https://en.wikipedia.org/wiki/ArcSDE" target="_blank">ArcSDE</a> or <a href="https://en.wikipedia.org/wiki/Oracle_Spatial_and_Graph" target="_blank">Oracle Spatial</a> at <a href="http://nvs-inc.com/" target="_blank">NVision</a> but over the last few years we have completely switched to PostGIS unless a client forces us to use a commercial solution.<br />
<br />
One of the main things I like about PostgreSQL is fantastic Python support with several options for database drivers. PostGIS is just a collection of PostgresSQL extensions and functions, so any database driver works. In this post, we're going to use the pure Python, dependency-free <a href="http://pg8000/" target="_blank">pg8000</a>, and <a href="https://github.com/GeospatialPython/pyshp" target="_blank">pyshp</a> modules to load a point shapefile into a posts database. Both modules can be quickly installed using <a href="https://github.com/pypa/pip" target="_blank">pip</a> or <a href="http://peak.telecommunity.com/DevCenter/EasyInstall" target="_blank">easy_install</a>.<br />
<br />
To get started, follow the PostGIS instructions to get up and running with a spatial database. In this example we'll use a database named "test" that is spatially enabled by PostGIS. We have pg8000 installed as our database interface and pushup to read a shapefile. The shapefile we'll use is a point shapefile containing information about some cities in the state of Mississippi. You can download the shapefile <a href="https://github.com/GeospatialPython/Learn/blob/master/MSCities_Geo_Pts.zip?raw=true" target="_blank">here</a>.<br />
<br />
The following code will set up a connection to our test geodatabase, create a table based on the fields of the shapefile, and then load the shapefile data into the table.<br />
<pre class="prettyprint"># Import the database driver and shapefile library
import pg8000
import shapefile
# Set up the datbase connection
connection = pg8000.connect(database="test", user="test",
password="test")
# Get the database cursor to execute queries
cursor = connection.cursor()
# Set up a shapefile reader for our shapefile
r = shapefile.Reader("MSCities_Geo_Pts/MSCities_Geo_Pts.shp")
# Build a query to create our "cities" table
# with an id and geometry column.
table_query = """CREATE TABLE cities (
id SERIAL,
PRIMARY KEY (id),
geom GEOMETRY(Point, 4326),
"""
# Get the shapefile fields but skip the first
# one which is a deletion flag used internally
# by dbf software
fields = r.fields[1:]
# We are going to keep track of the fields as a
# string. We'll use this query fragment for our
# insert query later.
field_string = ""
# Loop throug the fields and build our table
# creation query.
for i in range(len(fields)):
# get the field information
f = fields[i]
# get the field name and lowercase it for consistency
f_name = f[0].lower()
# add the name to the query and our string list of fields
table_query += f_name
field_string += f_name
# Get the proper type for each field. Note this check
# is not comprehensive but it convers the types in
# our sample shapefile.
if f[1] == "F":
table_query += " DOUBLE PRECISION"
elif f[1] == "N":
table_query += " INTEGER"
else:
table_query += " VARCHAR"
# If this field isn' the last, we'll add a comma
# and some formatting.
if i != len(fields)- 1:
table_query += ","
table_query += "\n"
table_query += " "
field_string += ","
# Close the query on the field.
else:
table_query += ")"
field_string += ",geom"
# Execute the table query
cursor.execute(table_query)
# Commit the change to the database or it won't stick.
# PostgreSQL is transactional which is good so nothing
# is stored until you're sure it works.
connection.commit()
# Create a python generator for the
# shapefiles shapes and records
shape_records = (shp_rec for shp_rec in r.iterShapeRecords())
# Loop through the shapefile data and add it to our new table.
for sr in shape_records:
# Get our point data.
shape = sr.shape
x, y = shape.points[0]
# Get the attribute data and set it up as
# a query fragment.
attributes = ""
for r in sr.record:
if type(r) == type("string"):
r = r.replace("'", "''")
attributes += "'{}'".format(r)
attributes += ","
# Build our insert query template for this shape record.
# Notice we are going to use a PostGIS function
# which can turn a WKT geometry into a PostGIS
# geometry.
point_query = """INSERT INTO cities
({})
VALUES ({}
ST_GEOMFROMTEXT('POINT({} {})',4326))"""
# Populate our query template with actual data.
format_point_query = point_query.format(field_string,
attributes, x, y)
# Insert the point data
cursor.execute(format_point_query)
# Everything went ok so let's update the database.
connection.commit()</pre>
<br />
The shapefile data is now available in PostGIS for loading into a GIS as a layer or performing spatial functions in PostGIS. Here's a screenshot of it loaded into QGIS as a PostGIS layer:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnhjEAWd76jFgrZjKw2WmHEWlZs4MJIYHlyRiee3Y4smEWX5v6CKTmaj1P0nFlU1Rn0I-13XTfL-R9SzPbCcaqb0iC9x_AAVtDIs0HMT3cNN7Nvr1B8LJprDB5yd66aSITEYjXe4w0hyphenhyphenY/s1600/Screen+Shot+2016-08-07+at+6.14.44+PM.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="250" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnhjEAWd76jFgrZjKw2WmHEWlZs4MJIYHlyRiee3Y4smEWX5v6CKTmaj1P0nFlU1Rn0I-13XTfL-R9SzPbCcaqb0iC9x_AAVtDIs0HMT3cNN7Nvr1B8LJprDB5yd66aSITEYjXe4w0hyphenhyphenY/s400/Screen+Shot+2016-08-07+at+6.14.44+PM.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
I've used this technique to download datasets from the web such as weather data, and update a table on an automated basis as new data is available. One thing we did not do in this example is add a spatial index. PostGIS can manage a lot of data - more than you can efficiently manage on the file system. Spatial indexing is a big part of that efficiency. It also allows multiple users to access that data at the same time. Once you load a data set, you typically do a query like the following to index the data:</div>
<pre class="prettyprint">index_query = """CREATE INDEX cities_geom_gix
ON cities USING GIST (geom);"""
</pre>
<br />
That's it! PostGIS has dozens of spatial functions that make analyzing data programmatically easy as well. With pg8000 and a PostgreSQL server, PostGIS becomes a giant pure Python GIS library with no compiled dependencies on your client platform.Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-63735732872698938122016-04-03T01:26:00.003-05:002016-08-22T09:58:23.216-05:00Pure Python NetCDF Reader<a href="https://thepythongischallenge.files.wordpress.com/2016/03/pynetcdf.png?w=990&h=486&crop=1" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="157" src="https://thepythongischallenge.files.wordpress.com/2016/03/pynetcdf.png?w=990&h=486&crop=1" width="320" /></a>Great news! Pure Python geospatial programmer Karim Bahgat created a depenency-free Python<b>pyncf</b>. If you're not familiar with NetCDF, it is technically a collection of software libraries and machine-independent data formats commonly used for multi-dimensional scientific data of any type. NetCDF currently uses a hierarchial data storage format called HDF5. NetCDF has a concept of dimensions, variables, and attributes and can scale from one-to-one relationships up to many-to-many relationships in all directions. In geospatial contexts, NetCDF is frequently used to store 2D and 3D raster or vector datasets over time. It is very popular for ocean observations and climatological data.<br />
NetCDF reader called<br />
<br />
He's posted the code and some more usage information on Github too: <a href="https://github.com/karimbahgat/pyncf">https://github.com/karimbahgat/pyncf</a><br />
<br />
Karim goes into detail on the ideas behind the library on his blog. But be careful, the examples on the blog aren't as current as the ones on Github: <a href="https://thepythongischallenge.wordpress.com/2016/03/26/pynetcdf-netcdf-files-in-pure-python/">https://thepythongischallenge.wordpress.com/2016/03/26/pynetcdf-netcdf-files-in-pure-python/</a><br />
<br />
As a quick test I installed pyncf using pip directly from Github:<br />
<span style="font-family: "courier new" , "courier" , monospace;">pip install https://github.com/karimbahgat/pyncf/archive/master.zip</span><br />
<br />
I then downloaded the sample time series orthogonal point data NetCDF sample file from here:<br />
<a href="https://www.nodc.noaa.gov/data/formats/netcdf/v1.1/">https://www.nodc.noaa.gov/data/formats/netcdf/v1.1/</a><br />
<br />
The actual link to the netcdf file is here:<br />
<a href="http://data.nodc.noaa.gov/thredds/fileServer/testdata/netCDFTemplateExamples/timeSeries/BodegaMarineLabBuoy.nc">/thredds/fileServer/testdata/netCDFTemplateExamples/timeSeries/BodegaMarineLabBuoy.nc</a><br />
<br />
What's nice about these NOAA samples is that they have a plain-text CDL (Common Data Language) version as well that let's you see what's in the file when you're experimenting with an API like this:<br />
<a href="http://data.nodc.noaa.gov/testdata/netCDFTemplateExamples/timeSeries/BodegaMarineLabBuoy.cdl">http://data.nodc.noaa.gov/testdata/netCDFTemplateExamples/timeSeries/BodegaMarineLabBuoy.cdl</a><br />
<br />
First, I imported the library and created a NetCDF object from the sample file:<br />
<br />
<pre class="prettyprint">import pyncf
nc = pyncf.NetCDF(filepath="BodegaMarineLabBuoy.nc")
</pre>
Then I created a variable to hold the header dictionary:<br />
<pre class="prettyprint">header = nc.header
</pre>
The header metadata contains a variety of summary data as nested dictionaries and list. We can access the product description by traversing that structure:<br />
<br />
<pre class="prettyprint">header["gatt_list"][1]["values"]
'These seawater data are collected by a moored fluorescence
and turbidity instrument operated at Cordell Bank, California, USA,
by CBNMS and BML. Beginning on 2008-04-23, fluorescence and turbidity
measurements were collected using a Wetlabs ECO Fluorescence
and Turbidity Sensor (ECO-FLNTUSB). The instrument depth of the
water quality sensors was 01.0 meter, in an overall water depth
of 85 meters (both relative to Mean Sea Level, MSL).
The measurements reflect a 10 minute sampling interval.'</pre>
<br />
The time values in this dataset are stored as epoch seconds. I accessed those and then converted the first and last into readable dates to see exactly what period this dataset spans:<br />
<br />
<pre class="prettyprint">import time
t = nc.read_dimension_values("time")
print time.ctime(t[0])
'Mon Jul 28 12:30:00 2008'
print time.time(t[-1])
'Wed Sep 10 10:31:00 2008'
</pre>
The pyncf codebase is considered an alpha version and is currently read only, but what a great addition to your pure Python geospatial toolbox! I wish this library was available when I updated "<a href="http://www.amazon.com/gp/product/1783552425/ref=as_li_tl?ie=UTF8&camp=1789&creative=9325&creativeASIN=1783552425&linkCode=as2&tag=geosppytho-20&linkId=HXO46G7K6XCZYOA3" target="_blank">Learning Geospatial Analysis with Python</a>"!<br />
<br />Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-20860548077414415672016-03-11T07:27:00.002-06:002016-03-11T07:27:39.580-06:00New Shapefile BrowserA developer named Adrián Eirís created a lightweight Python shapefile browser whose only compiled dependency is Tkinter which is included with most Python installations or is just a package manager command away. The tool is called "Ship-Shape-File-Navigator". It lets you preview a shapefile's geometry and attributes as well as perform common file operations like copying and renaming shapefiles as easily as if you were manipulating a single file (very handy!). So if you need to do some quick shapefile management and you don't wan to fire up QGIS or ArcMap, check it out. You can find out more at: <a href="http://shipshapefilenavigator.org/">http://ShipShapeFileNavigator.org</a><br />
<br />
<br />
Here's a quick screenshot showing the geometry preview:<br />
<div class="separator" style="clear: both; text-align: left;">
<a href="http://shipshp.bitbucket.org/images/screenshots/4_preview_preview.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://shipshp.bitbucket.org/images/screenshots/4_preview_preview.png" height="210" width="400" /></a></div>
<br />Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-82482523612647131322016-01-15T17:08:00.001-06:002016-01-15T17:08:20.107-06:002016 Book Sale!<div class="separator" style="clear: both; text-align: center;">
</div>
My publisher, Packt Publishing, started 2016 with a <span style="color: red;"><b>50% off sale</b></span> for the eBook versions of both the<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
first and second editions of "Learning Geospatial Analysis with Python"! Why would you want both? The first edition is all Python 2 and the second edition, besides additional code samples and other significant updates, is focused on Python 3. There are links on the right side of this page or you can just use code <b>LGAPS50 </b>at<a href="https://www.packtpub.com/" target="_blank"> https://www.packtpub.com</a>. This deal is good until Jan. 31 and I promise you won't be disappointed. <br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8e345i_8jraSK70423zDq1l4enTOvJeRvEyBAjoya0vqNcT40wvxZcCSI-Ve5kthSEVI1E-7jnY0JfoStnNqn1OzCqykTetBha-m3AtVMKcimG22P8PDgKYxLaZFDLQHhIFLLRERY32o/s1600/Snakes_04.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="150" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8e345i_8jraSK70423zDq1l4enTOvJeRvEyBAjoya0vqNcT40wvxZcCSI-Ve5kthSEVI1E-7jnY0JfoStnNqn1OzCqykTetBha-m3AtVMKcimG22P8PDgKYxLaZFDLQHhIFLLRERY32o/s200/Snakes_04.jpg" width="200" /></a></div>
<br />
When the first edition launched in 2013, there was one other geospatial python book out - Erik Westra's excellent "Geospatial Python Development". Now there are 12 on Packt's site alone and other publishers are working to catch up. I started this blog because I love Python and I love geospatial analysis but couldn't find much information out there. Now there is a wealth of excellent geospatial python info and software coming from all directions! <br />
<br />
You can help to continue to increase the use of Python in the geospatial field by discussing it online in blogs and forums, linking to interesting software and reading material on Twitter, releasing or contributing software on Github, and buying relevant books. Even in this digital age, books sales have a tremendous impact on the buzz and therefore the use of a given technology. And at the end of the day, mind share usually drives an industry to use a given technology. So however you participate, help spread the word that Python is the best language for geospatial technology (and pretty much everything else!). Of all the ways of "voting" for Python, buying books is the most point-and-click, dead simple, yet still personally rewarding way to make an impact!Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-31439945587757194302015-12-30T22:38:00.000-06:002016-02-10T15:14:42.153-06:00First Python 3 Geospatial Book<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: right; text-align: right;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs4sW-nW8ImpfGxr19UiKzho440INAbBEMAXk5PAPm8sfxsRcnv6GD47j14jdr2ZYS7C1xobiuh3eH8qv_MF_FvHLpXsSj9x1IwC5ngIlH9roKUL7X7hMFgdtFijatLcioSNunaIjvJPU/s1600/cover4.jpg" imageanchor="1" style="clear: right; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="208" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs4sW-nW8ImpfGxr19UiKzho440INAbBEMAXk5PAPm8sfxsRcnv6GD47j14jdr2ZYS7C1xobiuh3eH8qv_MF_FvHLpXsSj9x1IwC5ngIlH9roKUL7X7hMFgdtFijatLcioSNunaIjvJPU/s320/cover4.jpg" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><div style="text-align: left;">
<b>Chapter 8 demonstrates advanced geospatial modeling </b></div>
<div style="text-align: left;">
<b>using </b><b>Python 3 to create a color hillshade like t</b><b>his as well </b></div>
<div style="text-align: left;">
<b>as several other commonly used real-world examples.</b></div>
</td></tr>
</tbody></table>
The second edition of "Learning Geospatial Analysis with Python" is now available! To my<br />
<div style="text-align: right;">
</div>
knowledge, this is the first book dedicated to Python 3 and geospatial analysis. <br />
The book opens with the incredible story of how the grass-roots, open-source Ushahidi project stopped the Ebola epidemic in its tracks. It continues with one of the best histories of geospatial technologies you'll find. The chapters closes with you creating a complete GIS using only 60 lines of Python code!<br />
<br />
<b>Chapter 2</b> teaches you the data types you'll encounter on your journey into geospatial analysis. <b>Chapter 3 </b>surveys commonly-used geospatial software to give you the context you need to navigate geospatial solutions. <b>Chapter 4</b> covers Python geospatial tools. <b>Chapter 5</b> dives into GIS operations with Python. <b>Chapter 6 </b>teaches you Python 3 remote sensing techniques. And <b>Chapter 7</b> is a one-of-a-kind examination of elevation data with several examples. <b>Chapter 8</b> covers complete geospatial modeling techniques for real-world situations like the color hillshade in this post. <b>Chapter 9 </b>explores the emerging world of real-time tracking data. And <b>Chapter 10 </b>pulls the whole book together to build a complete application.<br />
<br />
You get over 100 examples with test data and code plus explanations in either a print or e-book available from Packt Publishing or pretty much any other online book seller. If you get the print book, or if you just want to see the images, I've made them available as a <a href="https://github.com/GeospatialPython/Learn/blob/master/LGAWP2ED_IMAGES.zip?raw=true" target="_blank">zip file download here</a>. <br />
<br />
While many books in this category do a great job of covering Python APIs to different libraries, this book focuses as much as possible on algorithms to show you the raw calculations that make up geospatial analysis. This approach allows you to not only run the examples but break things down line by line to see how it works and recombine pieces of code on your own to go beyond the book. <br />
<br />
Click the link on the right to check it out on Amazon!Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-34213451382782006822015-08-17T20:46:00.001-05:002015-08-17T20:56:53.031-05:00CSV to Shapefile<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWrk1Y5ET9YAZUDMQ2f3xkc_cMvMXhO9IVHG3sM1oO3SAddjjgIQz2K7yCBD1DtBYjIebwOGYJqrkGQZkpinwVx7trnP2q2Y_Np2WpBk94chJaStTdsi1LpDQTxFhfZNubzCXmV_ExHoo/s1600/csv2shp.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWrk1Y5ET9YAZUDMQ2f3xkc_cMvMXhO9IVHG3sM1oO3SAddjjgIQz2K7yCBD1DtBYjIebwOGYJqrkGQZkpinwVx7trnP2q2Y_Np2WpBk94chJaStTdsi1LpDQTxFhfZNubzCXmV_ExHoo/s320/csv2shp.png" width="308" /></a>Converting a <span style="font-family: inherit;">comma</span>-separated value (CSV) file to a shapefile comes up a lot on <a href="http://gis.stackexchange.com/search?q=csv+pyshp" target="_blank">forums</a>. This post is a quick example of transforming a CSV file to a shapefile using <a href="https://github.com/GeospatialPython/pyshp" target="_blank">PyShp</a> and the built-in Python <b>csv </b>module. CSV files usually contain fields separated by commas but might also use some other delimiter such as a pipe "|" or tab. The csv module defaults to commas but will let you set other delimiters.
Here is the file we'll be using that I slightly modified from a <a href="http://gis.stackexchange.com/questions/158560/converting-csv-files-into-shapefiles-using-python" target="_blank">GIS StackExchange</a> question flagged as duplicate:<br />
<br />
<pre><span style="font-family: Courier New, Courier, monospace;"><b>Name,Area,geometry
town,123,"POLYGON((73.953695297241 15.3555421329,73.951292037964 15.314816978733,73.986654281616 15.310015523128,73.982191085815 15.345775441645,73.953695297241 15.3555421329))"
forest,500,"POLYGON((73.938202857971 15.34362339739,73.944897651672 15.343375083164,73.943438529968 15.341554103151,73.942408561707 15.337084357624,73.941378593445 15.340395289421,73.937001228333 15.341471330955,73.938202857971 15.34362339739))"
lake,800,"POLYGON((73.995494842529 15.291305352455,73.999614715576 15.287165708411,74.000301361084 15.283357163686,73.997898101807 15.281701253096,73.99377822876 15.282363618901,73.99377822876 15.284516293317,73.992576599121 15.287496882943,73.992919921875 15.289815090018,73.995494842529 15.291305352455))"</b></span></pre>
<pre><span style="font-family: Courier New, Courier, monospace;">
</span></pre>
This notional CSV has three fields for the name of the feature, the area, and the geometry. The geometry consists of a WKT string describing the polygon coordinates. The polygons contain varying numbers of coordinates and the x,y values are separated by commas. None of this will confuse the csv module because they are contained in strings.
The code is very simple. We use the csv module to access the data as rows and fields. Then we manually parse the coordinates as strings. And finally we write the shapefile in a loop:
<br />
<pre><span style="font-family: Times, Times New Roman, serif;">
</span></pre>
<pre class="prettyprint">import csv
import shapefile
# Create a polygon shapefile writer
w = shapefile.Writer(shapefile.POLYGON)
# Add our fields
w.field("NAME", "C", "40")
w.field("AREA", "C", "40")
# Open the csv file and set up a reader
with open("sample.csv") as p:
reader = csv.DictReader(p)
for row in reader:
# Add records for each polygon for name and area
w.record(row["Name"], row["Area"])
# parse the coordinate string
wkt = row["geometry"][9:-2]
# break the coordinate string in to x,y values
coords = wkt.split(",")
# set up a list to contain the coordinates
part = []
# convert the x,y values to floats
for c in coords:
x,y = c.split(" ")
part.append([float(x),float(y)])
# create a polygon record with the list of coordinates.
w.poly(parts=[part])
# save the shapefile!
w.save("polys.shp")
</pre>
You can download the <a href="https://github.com/GeospatialPython/Learn/blob/master/csv2shp.py" target="_blank">code</a> and <a href="https://github.com/GeospatialPython/Learn/blob/master/sample.csv" target="_blank">CSV</a> file on GitHub.Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-3813585210878431432015-05-18T14:40:00.001-05:002015-05-18T15:01:07.906-05:00Basic Image Change Detection in QGIS and PythonChange detection allows you to automatically highlight the differences between two images in the same area if they are properly orthorectified. We'll do a simple difference change detection on two images, which are several years apart, to see the differences in urban development and the natural environment. Change detection algorithms can become quite sophisticated. However at the most basic level, you can do a simple, literal mathematical difference between the pixel values in the two images. That basic image difference is exactly what we'll do in this example. This recipe is from my book, the"<a href="http://www.amazon.com/gp/product/1783984988/ref=as_li_tl?ie=UTF8&camp=1789&creative=9325&creativeASIN=1783984988&linkCode=as2&tag=geosppytho-20&linkId=J4Y6EDB6KOFZQ26K" target="_blank">QGIS Python Programming Cookbook</a>". I show how do the same operation using the <b>gdalnumeric </b>module outside of QGIS in my other book "<a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20" target="_blank">Learning Geospatial Analysis with Python</a>".<br />
<br />
We're going to start with two aerial images taken several years apart. We'll subtract the before image from the after image and and visualize the result using a color ramp.<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg2Y9NxoJza1703KXxaKj9H_wpnQHXTWnq0YGFw0P6ri00J6c2N20RqSnYvv2k7l9rfjvcZb2BunQy9HvSp52qjCrQOjh3WJaxxpaptoHldB7i4obJBFPDcKqBQMLn7dt0hoAfZ9k9r1VU/s1600/change.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="161" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg2Y9NxoJza1703KXxaKj9H_wpnQHXTWnq0YGFw0P6ri00J6c2N20RqSnYvv2k7l9rfjvcZb2BunQy9HvSp52qjCrQOjh3WJaxxpaptoHldB7i4obJBFPDcKqBQMLn7dt0hoAfZ9k9r1VU/s400/change.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><b>The before and after aerial photos we'll use for this recipe. <br />Click for a larger version</b>.</td></tr>
</tbody></table>
<br />
<h3>
Getting ready</h3>
<br />
You will need a data directory to store the source images. I usually create a directory named <b>qgis_data</b> for these types of experiments. You can download the two images used in this example from <a href="https://github.com/GeospatialPython/qgis/blob/gh-pages/change-detection.zip?raw=true">https://github.com/GeospatialPython/qgis/blob/gh-pages/change-detection.zip?raw=true</a> and put them in a directory named change-detection in the rasters directory of your qgis_data directory. Note that the file is 55 megabytes, so it may take several minutes to download. Also, we'll step through this script a few lines at a time. It's meant to be run in the QGIS Python Console. You can download the entire script here: <a href="https://github.com/GeospatialPython/qgis/raw/master/4985_08_15-change.py">https://github.com/GeospatialPython/qgis/raw/master/4985_08_15-change.py</a><br />
<br />
<h3>
How to do it...</h3>
<br />
We'll use the QGIS raster calculator to subtract the images in order to get the difference, which will highlight significant changes. We'll also add a color ramp shader to the output in order to visualize the changes. To do this, we need to perform the following steps:<br />
<br />
First, we need to import the libraries that we need in to the QGIS console:<br />
<br />
<pre class="prettyprint">from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.analysis import *
</pre>
<br />
Now, we'll set up the path names and raster names for our images:<br />
<br />
<pre class="prettyprint">before = "/Users/joellawhead/qgis_data/rasters/change-detection/before.tif"
|after = "/Users/joellawhead/qgis_data/rasters/change-detection/after.tif"
beforeName = "Before"
afterName = "After"
</pre>
<br />
Next, we'll establish our images as raster layers:<br />
<br />
<pre class="prettyprint">beforeRaster = QgsRasterLayer(before, beforeName)
afterRaster = QgsRasterLayer(after, afterName)
</pre>
<br />
Then, we can build the calculator entries:<br />
<br />
<pre class="prettyprint">beforeEntry = QgsRasterCalculatorEntry()
afterEntry = QgsRasterCalculatorEntry()
beforeEntry.raster = beforeRaster
afterEntry.raster = afterRaster
beforeEntry.bandNumber = 1
afterEntry.bandNumber = 2
beforeEntry.ref = beforeName + "@1"
afterEntry.ref = afterName + "@2"
entries = [afterEntry, beforeEntry]
</pre>
<br />
Now, we'll set up the simple expression that does the math for remote sensing:<br />
<br />
<pre class="prettyprint">exp = "%s - %s" % (afterEntry.ref, beforeEntry.ref)</pre>
<br />
Then, we can set up the output file path, the raster extent, and pixel width and height:<br />
<br />
<pre class="prettyprint">output = "/Users/joellawhead/qgis_data/rasters/change-detection/change.tif"
e = beforeRaster.extent()
w = beforeRaster.width()
h = beforeRaster.height()
</pre>
<br />
Now, we perform the calculation:<br />
<br />
<pre class="prettyprint">change = QgsRasterCalculator(exp, output, "GTiff", e, w, h, entries)
change.processCalculation()
</pre>
<br />
Finally, we'll load the output as a layer, create the color ramp shader, apply it to the layer, and add it to the map, as shown here:<br />
<br />
<pre class="prettyprint">lyr = QgsRasterLayer(output, "Change")
algorithm = QgsContrastEnhancement.StretchToMinimumMaximum
limits = QgsRaster.ContrastEnhancementMinMax
lyr.setContrastEnhancement(algorithm, limits)
s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i = []
qri = QgsColorRampShader.ColorRampItem
i.append(qri(0, QColor(0,0,0,0), 'NODATA'))
i.append(qri(-101, QColor(123,50,148,255), 'Significant Itensity Decrease'))
i.append(qri(-42.2395, QColor(194,165,207,255), 'Minor Itensity Decrease'))
i.append(qri(16.649, QColor(247,247,247,0), 'No Change'))
i.append(qri(75.5375, QColor(166,219,160,255), 'Minor Itensity Increase'))
i.append(qri(135, QColor(0,136,55,255), 'Significant Itensity Increase'))
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)</pre>
<br />
You should end up with the following image. As a general rule of thumb, green shows area of change where something was built or trees were removed, and purple areas show where a road or structure was likely moved. There is lots of noise and false positives. But the algorithm does very well in areas of significant change.<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjc0xbHmE10lUR76d3pLS_NPMy1UuDBMOHSpY5OV-UwevYYVaZm1-dUmYGLBha97czMP_o1I5fSzzH6ApVi4Sg2C_lUJMGtSl60s8OYTIq2Ir8_9b3dM8zgat2h0FMy0e6WXk-kNCIOwwo/s1600/difference.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="291" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjc0xbHmE10lUR76d3pLS_NPMy1UuDBMOHSpY5OV-UwevYYVaZm1-dUmYGLBha97czMP_o1I5fSzzH6ApVi4Sg2C_lUJMGtSl60s8OYTIq2Ir8_9b3dM8zgat2h0FMy0e6WXk-kNCIOwwo/s400/difference.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><b>BAM! Change detection! Click for larger version</b></td></tr>
</tbody></table>
<h3>
How it works...</h3>
If a building is added in the new image, it will be brighter than its surroundings. If a building is removed, the new image will be darker in that area. The same holds true for vegetation, to some extent.<br />
<br />
Here is a close up of an area of intense change that demonstrates what we're talking about. In the following image, a single softball field was expanded into an entire softball field complex. In the upper left of the image you can see where the trees were cleared and bright-green grass and dirt infields were added. In the purple area in the lower right of the image you can see where the original field was rotated. Further in that direction, you can see some of the noise where the intensity changed between the images but it doesn't have much meaning. However if you zoomed out, the areas of significant change stand out very clearly which is usually the point of this kind of analysis:<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9QA70Y38lyvLx4PZZuZN03LgfHXAcsaUyas9bdOiSIw5HTFNY-oYFTaCN6tdCnkf052hCbxAm4NbXDKkJ6NhfeoxrNfrFnJoHuHXjVu74BgcgTCGlOtVKBOkY0JmILzFkHLewBQbde5k/s1600/COMPARISON.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="116" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9QA70Y38lyvLx4PZZuZN03LgfHXAcsaUyas9bdOiSIw5HTFNY-oYFTaCN6tdCnkf052hCbxAm4NbXDKkJ6NhfeoxrNfrFnJoHuHXjVu74BgcgTCGlOtVKBOkY0JmILzFkHLewBQbde5k/s400/COMPARISON.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><b>Looking at these softball fields close up really <br />shows what's happening. Click for larger version</b></td></tr>
</tbody></table>
<h3>
Summary</h3>
The concept is simple. We subtract the older image data from the new image data. Concentrating on urban areas tends to be highly-reflective and results in higher image pixel values.<br />
<br />Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-18833764543428880392015-05-13T10:44:00.001-05:002015-06-21T23:01:27.487-05:00Clipping a Shapefile in Pure PythonIn a previous post, I showed how to use pure Python to <a href="http://geospatialpython.com/2010/12/subsetting-shapefile-by-attributes.html" target="_blank">subset a shapefile by attributes</a>. That method is fast if the attributes define the data you want to clip. But sometimes you need to clip data spatially.<br />
<br />
In this post, we'll clip a road shapefile using bounding boxes. This is the simplest type of clip. In a future post, I'll demonstrate how to use irregular polygons which is more complex. However the ingredients I'll use for that recipe are already on this blog (this post + PyShp + the ray-casting point-in-polygon algorithm). <br />
<br />
Here's our goal for this example. We have a shapefile with the major roads of the United States which you can <a href="https://github.com/GeospatialPython/Learn/raw/master/roads.zip" target="_blank">download a sample dataset here</a>. We want to use a simple bounding box to subset the roads for Puerto Rico. Here's a view of part of that shapefile including Florida and Puerto Rico in the lower right of the image. It gives you a good idea of the density and number of records in this data.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuo_UA8PqZSGSNVY7P3KHNe9T3VEakErjnIXeZ5Q57cZLP67vVNXPmdMGeWky98C18kQWG2jor-oq1AjAKRlzQb9QMGmH3HFmvIsDoyf4b2c2FMPHb9WT5ZkNZOUR4iO1p4NkWTCY9PWE/s1600/roads.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuo_UA8PqZSGSNVY7P3KHNe9T3VEakErjnIXeZ5Q57cZLP67vVNXPmdMGeWky98C18kQWG2jor-oq1AjAKRlzQb9QMGmH3HFmvIsDoyf4b2c2FMPHb9WT5ZkNZOUR4iO1p4NkWTCY9PWE/s1600/roads.png" /></a></div>
<br />
<br />
<br />
Each road in the shapefile has its own shape record with the bounding box of the road in the record header. We'll just do a simple bounding box comparison to isolate the roads within our desired bounding box and then write those results to a new shapefile. You'll also need the <a href="https://raw.githubusercontent.com/GeospatialPython/pyshp/master/shapefile.py" target="_blank">latest version of PyShp</a> which provides the newer iterShapeRecords() method to provide and efficient iterator for both the shapes and the dbf attributes in a single object:<br />
<br />
<pre class="prettyprint">import shapefile
# Create a reader instance for our US Roads shapefile
r = shapefile.Reader("roadtrl020")
# Create a writer instance copying the reader's shapefile type
w = shapefile.Writer(r.shapeType)
# Copy the database fields to the writer
w.fields = list(r.fields)
# Our selection box that contains Puerto Rico
xmin = -67.5
xmax = -65.0
ymin = 17.8
ymax = 18.6
# Iterate through the shapes and attributes at the same time
for road in r.iterShapeRecords():
# Shape geometry
geom = road.shape
# Database attributes
rec = road.record
# Get the bounding box of the shape (a single road)
sxmin, symin, sxmax, symax = geom.bbox
# Compare it to our Puerto Rico bounding box.
# go to the next road as soon as a coordinate is outside the box
if sxmin < xmin: continue
elif sxmax > xmax: continue
elif symin < ymin: continue
elif symax > ymax: continue
# Road is inside our selection box.
# Add it to the new shapefile
w._shapes.append(geom)
w.records.append(rec)
# Save the new shapefile! (.shp, .shx, .dbf)
w.save("Puerto_Rico_Roads")
</pre>
<br />
And here's our result!<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhddZWW5buPWJq7EhMfGmjb8lQrAIx4bPXUVoEn0t5PVHkvT9eRYnUY4pldyXWHMffYm25Or1IFwV9emmMY7TISzep0Biycr3SnSCw0XlKj_CgZZBGX47hPywIjH3R6AX08weatyVhc6pg/s1600/PR_Roads.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhddZWW5buPWJq7EhMfGmjb8lQrAIx4bPXUVoEn0t5PVHkvT9eRYnUY4pldyXWHMffYm25Or1IFwV9emmMY7TISzep0Biycr3SnSCw0XlKj_CgZZBGX47hPywIjH3R6AX08weatyVhc6pg/s1600/PR_Roads.png" /></a></div>
Yes, there are dozens of tools and libraries that will do the same thing. But this is the only solution in the entire (known) universe that is implemented in pure Python. You don't need any compiled C code. If you have the ability to run Python 2 or 3, you can run this without any additional system permissions. It will run where most other libraries won't and when that's what your requirements are, nothing else will do!Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-89888694154383510462015-05-07T11:09:00.001-05:002015-05-07T11:11:14.225-05:00Geolocating Photos in QGIS<span style="font-family: FranklinGothic;">Photos taken with GPS-enabled cameras, including smartphones, store location information in
</span><span style="font-family: FranklinGothic;">the header of the file, in a format called EXIF tags. These tags are largely based on the same
header tags used by the TIFF image standard. If you have a directory full of photos, you can easily sort them by date or by name. But if you want to sort them by location, only a map will do. In this recipe, we'll use these tags to create
locations on a map in QGIS for some photos and provide links to open them. </span><br />
<div class="page" title="Page 40">
<div class="layoutArea">
<div class="column">
<h2>
<span style="color: #0b5394; font-family: FranklinGothic; font-size: large;">Getting ready</span></h2>
<span style="font-family: FranklinGothic;">You need to make sure you have <a href="http://www.qgis.org/en/site/">QGIS</a> installed as well as the <a href="http://effbot.org/zone/pil-index.htm">Python Imaging Library</a> (PIL) for the Python distribution that QGIS uses. QGIS should already have PIL included. You can download some sample geotagged images <a href="https://github.com/%20GeospatialPython/qgis/blob/gh-pages/photos.zip?raw=true">here</a> or you can use any photo taken outdoors with your smartphone. Unzip these photos and place them in a directory accessible to QGIS.</span><br />
<h2>
<span style="color: #0b5394; font-family: FranklinGothic; font-size: large;">How to do it...</span></h2>
PIL can parse EXIF tags. We will gather the file names of the photos, parse the location information, convert it to decimal degrees, create the point vector layer, add the photo locations, and add a QGIS action link to the attributes. To do this, we need to perform the following steps:<br />
<br />
In the QGIS Python Console, import the libraries that we'll need, including the Image module (PIL) and its ExifTags module for parsing image data and the glob module for doing wildcard file searches:<br />
<br />
<pre class="prettyprint">import glob
import Image
from ExifTags import TAGS
</pre>
<span style="font-family: FranklinGothic; font-size: x-small;"><br /></span>
<br />
Next, we'll create a function that can parse the header data:<br />
<br />
<pre class="prettyprint">def exif(img):
exif_data = {}
try:
i = Image.open(img)
tags = i._getexif()
for tag, value in tags.items():
decoded = TAGS.get(tag, tag)
exif_data[decoded] = value
except:
pass
return exif_data
</pre>
<br />
Now, we'll create a function that can convert degrees-minute-seconds to decimal degrees, which is how coordinates are stored in JPEG images:<br />
<br />
<pre class="prettyprint">def dms2dd(d, m, s, i):
sec = float((m * 60) + s)
dec = float(sec / 3600)
deg = float(d + dec)
if i.upper() == 'W':
deg = deg * -1
elif i.upper() == 'S':
deg = deg * -1
return float(deg)</pre>
<br />
Next, we'll define a function to parse the location data from the header data:
<br />
<br />
<pre class="prettyprint">def gps(exif):
lat = None
lon = None
if exif['GPSInfo']:
# Lat
coords = exif['GPSInfo']
i = coords[1]
d = coords[2][0][0]
m = coords[2][1][0]
s = coords[2][2][0]
lat = dms2dd(d, m ,s, i)
# Lon
i = coords[3]
d = coords[4][0][0]
m = coords[4][1][0]
s = coords[4][2][0]
lon = dms2dd(d, m ,s, i)
return lat, lon
</pre>
<br />
Next, we'll loop through the photos directory, get the file names, parse the location information, and build a simple dictionary to store the information, as follows:<br />
<br />
<pre class="prettyprint">photos = {}
photo_dir = "/Users/joellawhead/qgis_data/photos/"
files = glob.glob(photo_dir + "*.jpg")
for f in files:
e = exif(f)
lat, lon = gps(e)
photos[f] = [lon, lat]<span style="font-family: Times New Roman;"><span style="white-space: normal;">
</span></span></pre>
<br />
Now, we'll set up the vector layer for editing:<br />
<br />
<pre class="prettyprint">lyr_info = "Point?crs=epsg:4326&field=photo:string(75)"
vectorLyr = QgsVectorLayer(lyr_info, \"Geotagged Photos" , "memory")
vpr = vectorLyr.dataProvider()
</pre>
<br />
We'll add the photo details to the vector layer:<br />
<br />
<pre class="prettyprint">features = []
for pth, p in photos.items():
lon, lat = p
pnt = QgsGeometry.fromPoint(QgsPoint(lon,lat))
f = QgsFeature()
f.setGeometry(pnt)
f.setAttributes([pth])
features.append(f)
vpr.addFeatures(features)
vectorLyr.updateExtents()</pre>
<br />
Now, we can add the layer to the map and make the active layer:
<br />
<br />
<pre class="prettyprint">QgsMapLayerRegistry.instance().addMapLayer(vectorLyr)
iface.setActiveLayer(vectorLyr)
activeLyr = iface.activeLayer()
</pre>
<br />
Finally, we'll add an action that allows you to click on it and open the photo. Actions are a simple yet powerful feature of QGIS that let you define the result of certain user actions like clicking on features:
<br />
<br />
<pre class="prettyprint">actions = activeLyr.actions()
actions.addAction(QgsAction.OpenUrl, "Photos", \'[% "photo" %]')</pre>
<h2>
<span style="color: #0b5394; font-size: large;">
How it works...
</span></h2>
<div>
Using the included PIL EXIF parser, getting location information and adding it to a vector layer is relatively straight forward. This action is a default option for opening a URL. However, you can also use Python expressions as actions to perform a variety of tasks. The following screenshot shows an example of the data visualization and photo popup:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpCmRP84h199_jIf07io4-6rP5112xw-j6gMmRvAP030oXFjNJ_a4Rfu84mXbXP4QOiQbF2jCsZmXPoGrupwy0BdrKgeaM2huRJ1kp4AR0pIrO-hAIUNSjDj1umP-yTMEZABZ_wgQAOAw/s1600/4985_08_10.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="227" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpCmRP84h199_jIf07io4-6rP5112xw-j6gMmRvAP030oXFjNJ_a4Rfu84mXbXP4QOiQbF2jCsZmXPoGrupwy0BdrKgeaM2huRJ1kp4AR0pIrO-hAIUNSjDj1umP-yTMEZABZ_wgQAOAw/s400/4985_08_10.png" width="400" /></a></div>
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<h2>
<span style="color: #0b5394; font-size: large;">There's more...</span></h2>
</div>
<div>
You can download the complete code sample <a href="https://raw.githubusercontent.com/GeospatialPython/qgis/master/4985_08_14-geolocate.py">here</a>. Another plugin called Photo2Shape is available that performs a similar function, but it requires you to install an external EXIF tag parser while this approach uses the parser in PIL included with QGIS.</div>
</div>
</div>
</div>
Joel Lawheadhttp://www.blogger.com/profile/15645953215148786163noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-36683753881677497252014-12-19T14:19:00.001-06:002015-04-01T21:22:37.012-05:00Holiday SpecialPackt Publishing is offering a holiday special of $5 off all of their eBooks which include well over <br />
<div class="separator" style="clear: both; text-align: center;">
</div>
1,000 technical titles.<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://www.packtpub.com/sites/default/files/4985OS_2809_QGIS%20Python%20Programming.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="320" src="https://www.packtpub.com/sites/default/files/4985OS_2809_QGIS%20Python%20Programming.jpg" width="257" /></a></div>
<br />
<br />
<a href="https://www.packtpub.com/application-development/learning-geospatial-analysis-python">"Learning Geospatial Analysis with Python"</a> is <a href="https://www.packtpub.com/application-development/learning-geospatial-analysis-python">here</a>.<br />
<br />
Below are all of the titles at the moment that have geospatial-related content. More geospatial titles are on the way for 2104 including the little gem on the right :-)<br />
<br />
<ol>
<li><a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20" target="_blank">Learning Geospatial Analysis with Python</a></li>
<li><a href="http://www.amazon.com/gp/product/1849511543/ref=as_li_qf_sp_asin_il?ie=UTF8&tag=geosppytho-20&linkCode=as2&camp=1789&creative=9325&creativeASIN=1849511543" target="_blank">Python Geospatial Development</a></li>
<li><a href="http://www.amazon.com/dp/1849694443/ref=as_li_tf_til?tag=geosppytho-20&camp=0&creative=0&linkCode=as1&creativeASIN=1849694443&adid=1FAS3R0BF6XK25DF79Y6" target="_blank">Programming ArcGIS 10.1 Cookbook</a></li>
<li><a href="http://amzn.to/1gqc1rV" target="_blank">Applying and Extending Oracle Spatial</a></li>
<li><a href="http://amzn.to/1gC9ITP" target="_blank">Learning Bing Maps API</a></li>
<li><a href="http://amzn.to/OUJEaM" target="_blank">GeoServer Beginner’s Guide</a> </li>
<li><a href="http://amzn.to/1dvzJnh" target="_blank">Google Maps JavaScript API Cookbook</a></li>
<li><a href="http://amzn.to/1muYfI1" target="_blank">Python Data Visualization Cookbook</a></li>
<li><a href="http://amzn.to/1gCcMPz" target="_blank">OpenStreetMap</a></li>
<li><a href="http://amzn.to/1iAAme9" target="_blank">PostGIS Cookbook</a> </li>
<li><a href="http://amzn.to/OULhoK" target="_blank">OpenLayers 2.10 Beginner's Guide</a></li>
<li><a href="http://amzn.to/1maJd73" target="_blank">Learning QGIS 2.0</a></li>
<li><a href="https://www.packtpub.com/application-development/learning-qgis-second-edition">Learning QGIS - Second Edition</a></li>
<li><a href="https://www.packtpub.com/hardware-and-creative/geoserver-beginner%E2%80%99s-guide">GeoServer Beginner’s Guide</a></li>
<li><a href="https://www.packtpub.com/web-development/leafletjs-essentials">Leaflet.js Essentials</a></li>
</ol>
Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-54283781899395890202014-12-15T02:45:00.001-06:002014-12-19T11:37:56.489-06:00GPX2SHPHere's a nice Garmin GPX to shapefile converter using PyShp by Matt Rantala from his blog MapRantala.com:<br />
<div>
<br /></div>
<div>
<span style="-webkit-text-size-adjust: auto; background-color: rgba(255, 255, 255, 0);"><a href="http://www.maprantala.com/2014/05/16/garmin-gpx-to-shapefile-shp-conversion-gpx2shp-py/">http://www.maprantala.com/2014/05/16/garmin-gpx-to-shapefile-shp-conversion-gpx2shp-py/</a></span></div>
<div>
<br /></div>
<div>
<div class="separator" style="clear: both;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgFzbBwVgzJ-4Qy9jWYXRGrVIchszGhQ2-PuZMfYFned3DQ5itmwWsGYVtTw2aDJzq4gmZwbBlTZbYTUuTeQIDi2Xrq7F_jfuFca1NsOV399pNA84gkWibTxdKHqYASiaPxVzm2FZxEvi0/s640/blogger-image--577265626.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="179" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgFzbBwVgzJ-4Qy9jWYXRGrVIchszGhQ2-PuZMfYFned3DQ5itmwWsGYVtTw2aDJzq4gmZwbBlTZbYTUuTeQIDi2Xrq7F_jfuFca1NsOV399pNA84gkWibTxdKHqYASiaPxVzm2FZxEvi0/s320/blogger-image--577265626.jpg" width="320" /></a></div>
<div class="separator" style="clear: both;">
<b>Photo: FishStickTheatre.com</b></div>
<br /></div>
Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-71714274209180081702014-12-03T13:55:00.002-06:002014-12-19T14:03:19.573-06:00WKT EPSG Strings Made EasyThis quick post is a follow on to my posts on <a href="http://geospatialpython.com/2011/09/map-projections.html">map projections</a> and <a href="http://geospatialpython.com/2011/02/create-prj-projection-file-for.html">creating .prj files for shapefiles</a>.<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="http://spatialreference.org/static/img/sr_logo.jpg" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="http://spatialreference.org/static/img/sr_logo.jpg" height="53" width="200" /></a></div>
The .prj file format, created by Esri originally but now recognized by most geospatial software, is just a text file with a WKT string for the projection information. Most of projections can be referenced by a simple EPSG code which is just a number. Usually you'll have this number handy and need to look up the full WKT string on SpatialReference.org.<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0fxM2LCSKiK6hi92glVeI5TFONe7CyFRPhGuANAmj1jveajHwxDeWU6kRR2q3rvVHmMoCKi_uWQ0c0Iw71Pbv-6pqvB6SxPM_PHgZxfjrf3pZCxmmq1yxr5ZPssO1CyvuO43wqojIi78/s640/epsg-banner-1400-560-2.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0fxM2LCSKiK6hi92glVeI5TFONe7CyFRPhGuANAmj1jveajHwxDeWU6kRR2q3rvVHmMoCKi_uWQ0c0Iw71Pbv-6pqvB6SxPM_PHgZxfjrf3pZCxmmq1yxr5ZPssO1CyvuO43wqojIi78/s640/epsg-banner-1400-560-2.png" height="80" width="200" /></a></div>
<br />
<a href="http://www.vims.edu/people/forrest_dr/index.php">Dr. David Forrest</a>, a research scientist at the <a href="http://www.vims.edu/index.php">Virginia Institute of Marine Sciences</a>, sent me a simple function to make fetching WKT projection strings easy:<br />
<br />
<pre class="prettyprint">def getPRJwkt(epsg):
"""
Grab an WKT version of an EPSG code
usage getPRJwkt(4326)
This makes use of links like
http://spatialreference.org/ref/epsg/4326/prettywkt/
"""
import urllib
sr = "http://spatialreference.org"
f=urllib.urlopen(sr + "/ref/epsg/{0}/prettywkt/".format(epsg))
return (f.read())
</pre>
<br />
So then getting a WKT projection string is as easy as:<br />
<br />
<pre class="prettyprint">>>> wkt = getPRJwkt(4326)
>>> print wkt
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]
</pre>
<br />
There's another, newer website called <a href="http://epsg.io/">epsg.io</a> that has some additional features beyond spatialreference.org. The differences are in the FAQ on <a href="http://blog.klokantech.com/2014/03/epsgio-find-coordinate-systems-worldwide.html?utm_source=feedburner&utm_medium=feed&utm_campaign=Feed%3A+klokan-blog-osgeo+%28Blog%3A+Klokan+Petr+P%C5%99idal%29">this page</a>.Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-79678871688055363832014-12-01T23:22:00.000-06:002014-12-19T14:02:14.016-06:00ErrataIn Chapter 8 of my <a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20" target="">book</a>, the shapefile doesn't have a .prj file which isn't a big deal for the example<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjKCU0Jhpsp592HqO1hKcrDybHhxqZawztvgyXkJNny58NWAyaeFWPMJhz-jlki1EmaIlDHv0iBiRfVMijU1zYDVA3d3j2mFdc9kXoFrouHE8und6wWTHGa6BjcG_MbPvEPryXeYMhBTSk/s1600/1138_08_03.bmp" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjKCU0Jhpsp592HqO1hKcrDybHhxqZawztvgyXkJNny58NWAyaeFWPMJhz-jlki1EmaIlDHv0iBiRfVMijU1zYDVA3d3j2mFdc9kXoFrouHE8und6wWTHGa6BjcG_MbPvEPryXeYMhBTSk/s1600/1138_08_03.bmp" height="200" width="173" /></a></div>
but can be annoying if you try to open it in QGIS or ArcGIS. There was also an unused numpy array in the example so I updated the whole package. You can download the updated script and data here:<br />
<a href="https://geospatialpython.googlecode.com/svn/NDVI-update.zip">https://geospatialpython.googlecode.com/svn/NDVI-update.zip</a><br />
<br />
This example is the very first example in the chapter. The second example shows you how to classify the output data to create the image shown here.Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0tag:blogger.com,1999:blog-2044555433406398156.post-51333136010898230552014-12-01T11:37:00.004-06:002014-12-01T11:37:35.215-06:00"Learning Geospatial Analysis with Python" Cyber Monday SaleThe print edition of <b><a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20">Learning Geospatial Analysis with Python</a></b> qualifies for a 30% discount on Amazon.com until midnight tonight using discount code <span style="font-family: Courier New, Courier, monospace;"><span style="color: red;"><b>HOLIDAY30</b></span> </span>at checkout. This is for the print edition only which is black and white as are most programming books as opposed to the eBook which is full color. Very soon I will make an image catalog available for download for people with the print edition. This detail isn't critical as most of the illustrations are not color dependent with the exception of a 1 or 2 remote sensing examples where color is helpful. The image catalog will give you the best of both worlds and with a 30% savings to boot!<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuJfUxXcm794JfC8MAwJbwNfjNGTPphAhqgx8gAHIySi9sD-0ls1iy4s-LDKoaRJH68aasv3iM-BsoDgOGCXHK_VpDBkx5MDnPQsaXt0Ff83xZw7uBANzpVYAUDoT8dZ55WhC4b9IKuMA/s1600/1138OS_Cov.jpg" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"></a><a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20"><img border="0" src="http://nicktumminello.com/wp-content/uploads/2013/12/30-percent-off-discount-sale-icon_2.png" height="320" width="262" /></a> <a href="http://www.amazon.com/gp/product/1783281138/ref=as_li_ss_tl?ie=UTF8&camp=1789&creative=390957&creativeASIN=1783281138&linkCode=as2&tag=geosppytho-20"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuJfUxXcm794JfC8MAwJbwNfjNGTPphAhqgx8gAHIySi9sD-0ls1iy4s-LDKoaRJH68aasv3iM-BsoDgOGCXHK_VpDBkx5MDnPQsaXt0Ff83xZw7uBANzpVYAUDoT8dZ55WhC4b9IKuMA/s1600/1138OS_Cov.jpg" height="320" width="259" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<table><tbody>
<tr>
<td><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEAnK_D6_TIfbNGu0yKK_Aurs9SOoysyAfH26tNYMKi6uNk_xThj9mPiRbYtzRtW6WWchA2vWsXaQkiNMlGAeNZdCqk4s-kPCHa5_SX6Rq7z_WXrQxbnkOvXbSW78dUPIsEUmfULcZsCY/s1600/1138_07_04.jpg" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEAnK_D6_TIfbNGu0yKK_Aurs9SOoysyAfH26tNYMKi6uNk_xThj9mPiRbYtzRtW6WWchA2vWsXaQkiNMlGAeNZdCqk4s-kPCHa5_SX6Rq7z_WXrQxbnkOvXbSW78dUPIsEUmfULcZsCY/s1600/1138_07_04.jpg" height="157" width="200" /></a>
</td>
<td><a href="http://1.bp.blogspot.com/-tBEZzcVWpjU/Up9tD_w5rAI/AAAAAAAAAP4/T6IadXfPZmg/s1600/1138_07_01.jpg" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="http://1.bp.blogspot.com/-tBEZzcVWpjU/Up9tD_w5rAI/AAAAAAAAAP4/T6IadXfPZmg/s1600/1138_07_01.jpg" height="158" width="200" /></a>
</td>
</tr>
</tbody></table>
Joel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.com0