1 / 19

Raster

Python Geoprocessing Outside the ESRI Box. Raster. GDAL. Derived data. Vector. Python. OGR. -Microsoft ADO -others. Relational DB. Why Use non-ESRI libraries?. Large performance gains for some operations – more scalable Greater flexibility, if willing to program Full customization

menora
Download Presentation

Raster

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. Python Geoprocessing Outside the ESRI Box Raster GDAL Deriveddata Vector Python OGR -Microsoft ADO -others Relational DB

  2. Why Use non-ESRI libraries? • Large performance gains for some operations – more scalable • Greater flexibility, if willing to program • Full customization • Large number formats/data sources supported • Open source

  3. GDAL – Geospatial Data Abstraction Library • www.gdal.org • Used in several major software packages • Installation – FWtools Windows binary • Includes a Python instance • Many formats supported: • Arc GRID, SDE raster, Erdas Imagine, GeoTiff, NetCDF, GRASS, several others • Command-line utilities:gdalinfo, gdal_translate, gdal_merge.py, gdal_rasterize • Call from Python - os.system(command)

  4. Programming with GDAL • Library presents a single abstract data model for supported formats • Classes include: • GDALDriver • GDALDataset • GDALRasterBand • others

  5. Open a Raster Dataset for Reading import gdal from gdalconst import * # open GDAL dataset gdata = gdal.Open(‘e:\\data\\z02_evt.img’, GA_ReadOnly) if gdata is None: print 'could not open the raster' sys.exit(1)

  6. Get Projection, Extent, and Cell Size # projection in OpenGIS WKT format proj = gdata.GetProjection() # extent, cell size xSize = gdata.RasterXSize ySize = gdata.RasterYSize geotr = gdata.GetGeoTransform() if not geotr is None: originX = float(geotr[0]) originY = float(geotr[3]) pixelXSize = abs(float(geotr[1])) pixelYSize = abs(float(geotr[5])) geotr = None

  7. Read the First Row of Pixels into a Numeric Array import Numeric # fetch the first band of the raster band = gdata.GetRasterBand(1) # upper left pixel is 0,0 # read first row into an array row = band.ReadAsArray(0, 0, xSize, 1)

  8. Overlay Plots on 30-m Imagery and Other Predictors

  9. Overlay Plots on Vegetation Layers for Accuracy Assessment

  10. Point Overlay Application • Batches of 120 raster layers Imagine (.img) and Arc GRID • 1,000s to 10,000s points per batch • Several batches per week • Need to support single-pixel, as well as multiple-pixel kernels (e.g., 3x3, 5x5, or irregular kernels) • Need access to the raw pixel values • Need to process batches in minutes, not hours, to support production environment and rapid experimentation

  11. OGR Simple Feature Library • Low-level access to vector data • OGRSpatialReference • OGRCoordinateTransformation • Utility programs include ogr2ogr: • Transform formats • Can write kml

  12. Microsoft ADO • Single abstract data model for nearly all database formats • Proprietary, but built in on Windows • Powerful and easy to use • COM objects – instantiate like GP • mainly Connection, Recordset

  13. Microsoft ADO import win32com.client # Create Connection, Recordset objects cn = win32com.client.Dispatch(r"ADODB .Connection") rs = win32com.client.Dispatch(r"ADODB .Recordset") # Open connection to Access database cn.Open(r'Provider=Microsoft.Jet.OLEDB. 4.0;Data Source=' + path_to_mdb)

  14. Microsoft ADO # Execute SQL directly with Connection # sql = an SQL statement (“SELECT …”) cn.Execute(sql) # Or, open a Recordset rs.Open(sql, cn, 3, 1) # Loop over all records rs.MoveFirst() while not rs.EOF: # do stuff rs.MoveNext() rs.Close()

  15. More resources… ArcUser, April-June 2005 http://www.esri.com/news/arcuser/0405/files/python.pdf

  16. More resources…

  17. More resources… Python and ADO: http://www.markcarter.me.uk/computing/python/ado.html Python Cookbook: http://aspn.activestate.com/ASPN/Python/Cookbook/

  18. Chris Toney RMRS Missoula Fire Sciences Lab christoney@fs.fed.us

  19. >>>import this

More Related