| Hosted by CoCalc | Download

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))