1 / 29

Introduction to Open Source RS/GIS programming with Python

Introduction to Open Source RS/GIS programming with Python. Chris Garrard RS/GIS Laboratory Utah State University. Why bother with FOSS?. Free and Open Source Software (FOSS) Affordable (free, as in free beer) Fits the budget of small companies or individuals Open (as in free speech)

annice
Download Presentation

Introduction to Open Source RS/GIS programming with Python

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Introduction to Open Source RS/GIS programming with Python Chris Garrard RS/GIS Laboratory Utah State University

  2. Why bother with FOSS? Free and Open Source Software (FOSS) • Affordable (free, as in free beer) • Fits the budget of small companies or individuals • Open (as in free speech) • Use it for whatever you want • Distribute it if you want • Modify it if you want

  3. Why bother with FOSS? • Helpful user community • Including the developers themselves • Fast bug fixes • Not limited to Windows • Lightweight

  4. Cons • Lightweight • Not the ESRI geoprocessor • Smaller user community than some other packages

  5. Free != inferior • Apache, MySQL, Linux, Firefox, etc. • ESRI likes Python • GDAL is used by popular commercial packages • ArcGIS, Erdas Imagine, ER Mapper

  6. Python • Free, even for resale • Powerful and has extensive libraries • Easy to integrate with other tools • Runs on many different operating systems • Easy to learn

  7. GDAL and OGR libraries • Free, even for resale • Available for many different operating systems • Written in ANSI C and C++ • Bindings exist for Python, C#, Perl, Java, Ruby, R, VB6

  8. OGR • OGR Simple Features Library • Vector data access • Software-specific formats • shapefiles, personal geodatabases, ArcSDE, MapInfo, GRASS, Microstation • Open formats • TIGER/Line, SDTS, GML, KML • Databases • MySQL, PostgreSQL, Oracle Spatial, Informix, ODBC

  9. GDAL • Geospatial Data Abstraction Library • Raster data access • Supports about 100 different formats • ArcInfo grids, ArcSDE raster, Imagine, Idrisi, ENVI, GRASS, GeoTIFF • HDF4, HDF5 • USGS DOQ, USGS DEM • ECW, MrSID • TIFF, JPEG, JPEG2000, PNG, GIF, BMP

  10. Real-world examples • Drill through a few hundred raster images and pull out cell values at several hundred points (~3.5 million values in output) • Batch compute NDWI (Normalized Difference Water Index) for Landsat imagery • Batch clip extraneous data from edges of images

  11. Real-world examples • Created a text file with info (bounding coordinates, dimensions, etc) about several hundred MSS images which was then used to input data into a database • Created polygons showing bounding coordinates of actual data in several hundred raster images • Logistic regression

  12. Location Calculator

  13. <%@ Language=Python %> def getdata (keys=''): ''' possible ways to call this function: value=getdata('key') dict=getdata(('key1','key2')) #get subset dict=getdata() #return everything It assumes you don't have the same key for GET and POST methods ''' import types key_type=type(keys) d_data={} #initialize dictionary if keys=='': #if they don't supply keys then return everything for key in Request.Form: d_data[key]=Request.Form(key) for key in Request.QueryString: d_data[key]=Request.QueryString(key) return d_data elif key_type == types.StringType: #if they provide a single string value=Request.Form(keys) if (value != 'None') and (str(value) == 'None'): return Request.QueryString(keys) else: return value #if they provide a list then return a dictionary with all the key/values elif key_type == types.TupleType or key_type == types.ListType: for key in keys: value=Request.Form(key) #now check if the data was empty, if so look at QueryString if (value != 'None') and (str(value) == 'None'): value=Request.QueryString(key) data[key]=value return d_data if str(key) == 'utm': X_UTM=getdata('UTMx') Y_UTM=getdata('UTMy') datum=getdata('datumutm')

  14. UTM NAD27 Zone 12 North 420949.17 4513683.50 Projects to:

  15. UTM NAD27 Zone 12 North 420949.17 4513683.50 Projects to:

  16. UTM NAD27 Zone 12 North 420949.17 4513683.50 Projects to:

  17. Geocoding

  18. Shapefile from SQL

  19. Ex. 1: Reading vector data • Open a point shapefile • For each point: • Get the “id” and “cover” attributes • Get the x,y coordinates for the point • Print “id x y cover”

  20. Ex. 2: Extracting vector data • Extract all “Quacking aspen woodland alliance” SWReGAP field sites in Wasatch county from the statewide data set • Open a point shapefile with SWReGAP sample sites • Open a Utah counties polygon shapefile • Create a new point shapefile • Use the county file to limit the SWReGAP shapefile to Wasatch County

  21. Ex. 2: Extracting vector data • Limit SWReGAP to “Quaking aspen woodland alliance” • For each point in SWReGAP (remember it is now filtered spatially and by attribute): • Copy the feature into the new shapefile • Close all of the files

  22. Ex. 3: Generate Near Table tool • Given a point shapefile (SWReGAP field sites) and a line shapefile (roads), calculate the distance from each point to the nearest road • Open both shapefiles • For each point: • Limit the search to the provided radius by buffering the point and using that to set a spatial filter on the roads

  23. Ex. 3: Generate Near Table tool • For each point, continued: • Get the distance to the first line and its FID • Loop through the rest of the lines and if the distance is less than the “remembered” distance, then remember the new distance and FID • Write out the point FID and the remembered line FID and distance to a text file • Close the shapefiles and the text file • Needs an ArcInfo license in ArcMap (can probably write your own script, though)

  24. Ex. 4: Extract pixel values • Open a point shapefile • Open a 3-band raster file • For each point in the shapefile: • Print out the pixel values at that location for all three bands • “Sample” tool requires Spatial Analyst (to be fair, I think I may have done this with the arcgisscripting module without SA)

  25. Ex. 5: Compute NDVI • Create a Normalized Difference Vegetation Index for an ASTER image • Open the image • Get the image dimensions and block size • Use this info to create a new image • Read in one block at a time from the image • Compute the NDVI: (band3-band2)/(band3+band2) • Write the block to the new image • Compute the statistics for the new image

  26. Ex. 5: Compute NDVI • Set the georeferencing and projection for the new image • Build pyramids for the new image • Needs the Spatial Analyst extension, and I don’t know if you could do it with anything in the geoprocessor without SA – I’m sure you could do it with ArcObjects

  27. Ex. 6: NDVI in ArcMap • Run the NDVI script in ArcMap • Change the script to accept parameters • Change the script to send error messages to the geoprocessor • Add the script to a toolbox

  28. Thanks! • Course materials at www.gis.usu.edu/~chrisg/python • I can be contacted at chrisg@gis.usu.edu

More Related