Sharedsupport / 2014-11-30-gdal.sagewsOpen in CoCalc
Authors: Harald Schilly, ℏal Snyder, William A. Stein
License: GNU General Public License v3.0
Description: Examples for support purposes.

GDAL: Geospatial Data Abstraction Library

See http://www.gdal.org/

Documentation for using GDAL from Python:

http://pcjericks.github.io/py-gdalogr-cookbook/
# If you get an error when you evaluate this, restart your project server in project settings. import gdal from osgeo import ogr, osr, gdal
import sys from osgeo import gdal version_num = int(gdal.VersionInfo('VERSION_NUM')) version_num
1110100
from osgeo import gdal # Enable GDAL/OGR exceptions gdal.UseExceptions() # open dataset that does not exist ds = gdal.Open('test.tif') # results in python RuntimeError exception that # `test.tif' does not exist in the file system
Error in lines 5-7 Traceback (most recent call last): File "/projects/4a5f0542-5873-4eed-a85c-a18c706e8bcd/.sagemathcloud/sage_server.py", line 865, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> RuntimeError: `test.tif' does not exist in the file system, and is not recognised as a supported dataset name.
# IMPORTANT! In Sage you must do this to avoid type issues. RealNumber = float; Integer = int
# Create a point from osgeo import ogr point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(1198054.34, 648493.09) print point.ExportToWkt()
POINT (1198054.34 648493.09 0)
# Create a LineString from osgeo import ogr line = ogr.Geometry(ogr.wkbLineString) line.AddPoint(1116651.439379124, 637392.6969887456) line.AddPoint(1188804.0108498496, 652655.7409537067) line.AddPoint(1226730.3625203592, 634155.0816022386) line.AddPoint(1281307.30760719, 636467.6640211721) print line.ExportToWkt()
LINESTRING (1116651.439379123970866 637392.696988745592535 0,1188804.010849849553779 652655.740953706670552 0,1226730.362520359223709 634155.08160223858431 0,1281307.307607189984992 636467.66402117209509 0)
from osgeo import ogr # Create ring ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(1179091.1646903288, 712782.8838459781) ring.AddPoint(1161053.0218226474, 667456.2684348812) ring.AddPoint(1214704.933941905, 641092.8288590391) ring.AddPoint(1228580.428455506, 682719.3123998424) ring.AddPoint(1218405.0658121984, 721108.1805541387) ring.AddPoint(1179091.1646903288, 712782.8838459781) # Create polygon poly = ogr.Geometry(ogr.wkbPolygon) poly.AddGeometry(ring) print poly.ExportToWkt()
0 POLYGON ((1179091.164690328761935 712782.883845978067257 0,1161053.021822647424415 667456.268434881232679 0,1214704.933941904921085 641092.828859039116651 0,1228580.428455505985767 682719.312399842427112 0,1218405.065812198445201 721108.180554138729349 0,1179091.164690328761935 712782.883845978067257 0))
# Polygon with holes from osgeo import ogr # Create outer ring outRing = ogr.Geometry(ogr.wkbLinearRing) outRing.AddPoint(1154115.274565847, 686419.4442701361) outRing.AddPoint(1154115.274565847, 653118.2574374934) outRing.AddPoint(1165678.1866605144, 653118.2574374934) outRing.AddPoint(1165678.1866605144, 686419.4442701361) outRing.AddPoint(1154115.274565847, 686419.4442701361) # Create inner ring innerRing = ogr.Geometry(ogr.wkbLinearRing) innerRing.AddPoint(1149490.1097279799, 691044.6091080031) innerRing.AddPoint(1149490.1097279799, 648030.5761158396) innerRing.AddPoint(1191579.1097525698, 648030.5761158396) innerRing.AddPoint(1191579.1097525698, 691044.6091080031) innerRing.AddPoint(1149490.1097279799, 691044.6091080031) # Create polygon poly = ogr.Geometry(ogr.wkbPolygon) poly.AddGeometry(outRing) poly.AddGeometry(innerRing) print poly.ExportToWkt()
0 0 POLYGON ((1154115.274565846892074 686419.444270136067644 0,1154115.274565846892074 653118.257437493419275 0,1165678.186660514445975 653118.257437493419275 0,1165678.186660514445975 686419.444270136067644 0,1154115.274565846892074 686419.444270136067644 0),(1149490.109727979870513 691044.609108003089204 0,1149490.109727979870513 648030.576115839648992 0,1191579.109752569813281 648030.576115839648992 0,1191579.109752569813281 691044.609108003089204 0,1149490.109727979870513 691044.609108003089204 0))
# Multipoint from osgeo import ogr multipoint = ogr.Geometry(ogr.wkbMultiPoint) point1 = ogr.Geometry(ogr.wkbPoint) point1.AddPoint(1251243.7361610543, 598078.7958668759) multipoint.AddGeometry(point1) point2 = ogr.Geometry(ogr.wkbPoint) point2.AddPoint(1240605.8570339603, 601778.9277371694) multipoint.AddGeometry(point2) point3 = ogr.Geometry(ogr.wkbPoint) point3.AddPoint(1250318.7031934808, 606404.0925750365) multipoint.AddGeometry(point3) print multipoint.ExportToWkt()
0 0 0 MULTIPOINT (1251243.736161054344848 598078.795866875909269 0,1240605.857033960288391 601778.927737169433385 0,1250318.703193480847403 606404.092575036454946 0)
from osgeo import ogr multiline = ogr.Geometry(ogr.wkbMultiLineString) line1 = ogr.Geometry(ogr.wkbLineString) line1.AddPoint(1214242.4174581182, 617041.9717021306) line1.AddPoint(1234593.142744733, 629529.9167643716) multiline.AddGeometry(line1) line1 = ogr.Geometry(ogr.wkbLineString) line1.AddPoint(1184641.3624957693, 626754.8178616514) line1.AddPoint(1219792.6152635587, 606866.6090588232) multiline.AddGeometry(line1) print multiline.ExportToWkt()
0 0 MULTILINESTRING ((1214242.417458118172362 617041.971702130627818 0,1234593.142744733020663 629529.916764371562749 0),(1184641.36249576928094 626754.817861651419662 0,1219792.615263558691368 606866.609058823203668 0))
from osgeo import ogr multipolygon = ogr.Geometry(ogr.wkbMultiPolygon) # Create ring #1 ring1 = ogr.Geometry(ogr.wkbLinearRing) ring1.AddPoint(1204067.0548148106, 634617.5980860253) ring1.AddPoint(1204067.0548148106, 620742.1035724243) ring1.AddPoint(1215167.4504256917, 620742.1035724243) ring1.AddPoint(1215167.4504256917, 634617.5980860253) ring1.AddPoint(1204067.0548148106, 634617.5980860253) # Create polygon #1 poly1 = ogr.Geometry(ogr.wkbPolygon) poly1.AddGeometry(ring1) multipolygon.AddGeometry(poly1) # Create ring #2 ring2 = ogr.Geometry(ogr.wkbLinearRing) ring2.AddPoint(1179553.6811741155, 647105.5431482664) ring2.AddPoint(1179553.6811741155, 626292.3013778647) ring2.AddPoint(1194354.20865529, 626292.3013778647) ring2.AddPoint(1194354.20865529, 647105.5431482664) ring2.AddPoint(1179553.6811741155, 647105.5431482664) # Create polygon #2 poly2 = ogr.Geometry(ogr.wkbPolygon) poly2.AddGeometry(ring2) multipolygon.AddGeometry(poly2) print multipolygon.ExportToWkt()
0 0 0 0 MULTIPOLYGON (((1204067.054814810631797 634617.598086025333032 0,1204067.054814810631797 620742.10357242426835 0,1215167.450425691669807 620742.10357242426835 0,1215167.450425691669807 634617.598086025333032 0,1204067.054814810631797 634617.598086025333032 0)),((1179553.681174115510657 647105.543148266384378 0,1179553.681174115510657 626292.30137786467094 0,1194354.208655290072784 626292.30137786467094 0,1194354.208655290072784 647105.543148266384378 0,1179553.681174115510657 647105.543148266384378 0)))
import ogr # Given a test polygon poly_Wkt= "POLYGON((-107.42631019589980212 40.11971708125970082,-107.42455436683293613 40.12061219666851741,-107.42020981542387403 40.12004414402532859,-107.41789122063043749 40.12149008687303819,-107.41419947746419439 40.11811617239460048,-107.41915181585792993 40.11761695654455906,-107.41998470913324581 40.11894245264452508,-107.42203317637793702 40.1184088144647788,-107.42430674991324224 40.1174448122981957,-107.42430674991324224 40.1174448122981957,-107.42631019589980212 40.11971708125970082))" geom_poly = ogr.CreateGeometryFromWkt(poly_Wkt) print geom_poly
POLYGON ((-107.426310195899802 40.119717081259701,-107.424554366832936 40.120612196668517,-107.420209815423874 40.120044144025329,-107.417891220630437 40.121490086873038,-107.414199477464194 40.1181161723946,-107.41915181585793 40.117616956544559,-107.419984709133246 40.118942452644525,-107.422033176377937 40.118408814464779,-107.424306749913242 40.117444812298196,-107.424306749913242 40.117444812298196,-107.426310195899802 40.119717081259701))