python gdal tutorial: geometry and projection

Create an empty geometry object: ogr.Geometry

The methods used to define various geometries are different (point, line, polygon, etc)

To create a new point, use the method AddPoint( , < ;y>, []). The z coordinate is generally omitted, and the default value is 0

For example:

point = ogr.Geometry(ogr.wkbPoint)

point.AddPoint(10,20)

New line

Use AddPoint(< x>, , []) to add a point

Use SetPoint(, , , []) to change the coordinates of a point

For example, the following code changes the coordinates of point 0:

line = ogr.Geometry(ogr.wkbLineString)

line.AddPoint(10,10)

line.AddPoint(20,20)

line. SetPoint(0,30,30) #(10,10) -> (30,30)

Count the number of all points

print line.GetPointCount()

Read the x coordinate and y coordinate of point 0

print line.GetX(0)

print line.GetY(0)

To create a new polygon, first create a new ring, and then add the ring to the polygon object.

How to create a ring? First create a new ring object, and then add points to it one by one.

ring = ogr.Geometry(ogr.wkbLinearRing)

ring.AddPoint(0,0)

ring.AddPoint(100,0)

ring.AddPoint(100,100)

ring.AddPoint(0,100)

At the end, use CloseRings to close the ring, or set the coordinates of the last point to be the same as the first point.

ring.CloseRings()

ring.AddPoint(0,0)

The following is an example to create a box. This is a polygon object, composed of two layers of rings. Outring = OGR.Geometry (OGR.WKBLINEARING)

Outring.addpoint (0,0). ) outring.AddPoint(0,0)

inring = ogr.Geometry(ogr.wkbLinearRing)inring = ogr.Geometry(ogr.wkbLinearRing)

inring.AddPoint(25,25)

inring.AddPoint(75,25)

inring.AddPoint(75,75)

inring.AddPoint(25,75)

inring.CloseRings()

polygon = ogr.Geometry(ogr.wkbPolygon)

polygon.AddGeometry(outring)

polygon .AddGeometry(inring)

The last three sentences are more important, that is, first create a polygon object, and then add the outer ring and inner ring

The following sentence can help you count how many rings your polygon can have

print polygon.GetGeometryCount()

Read ring from polygon, the order of index is the same as the order of adding ring when creating polygon

outring = polygon.GetGeometryRef(0)

inring = polygon.GetGeometryRef(1)

Create multi geometry

such as MultiPoint, MultiLineString, MultiPolygon

Use AddGeometry to add ordinary geometric shapes to the composite geometry, for example:

multipoint = ogr.Geometry(ogr.wkbMultiPoint)

point = ogr. Geometry(ogr.wkbPoint)point = ogr.Geometry(ogr.wkbPoint)

point.AddPoint(10,10)

multipoint.AddGeometry(point)

point.AddPoint(20,20)

multipoint.AddGeometry( point)

Reading Geometry in MultiGeometry is the same as reading ring from Polygon. It can be said that Polygon is a built-in MultiGeometry.

Don’t delete an existing Feature’s Geometry, it will crash python.

You can only delete the Geometry created during the running of the script, for example, manually created, or automatically created by calling other functions. Even if this Geometry has been used to create other Features, you can still delete it.

For example: Polygon.Destroy()

Regarding projection Projections, use SpatialReference object

A variety of Projections, GDAL supports WKT, PROJ.4, ESPG, USGS, ESRI.prj

can be read from layer and Geometry Get Projections, for example:

spatialRef = layer.GetSpatialRef()

spatialRef = geom.GetSpatialReference()

Projection information is generally stored in the .prj file. If there is no such file, the above function returns None

Create a new one Projection:

First import the osr library, then use osr.SpatialReference() to create a SpatialReference object

Then use the following statements to import the projection information to the SpatialReference object

  •ImportFromWkt()

  •ImportFromEPSG( ; ;)

 •ImportFromUSGS(< proj_code>, )

•ImportFromXML()

Export Projection, use the following statement to export it as a string

​ •ExportToWkt()

​ •ExportToPrettyWkt()

​ •ExportToPro​ j4()

  •ExportToPCI()

  •ExportToUSGS()

•ExportToXML()

To perform projection transformation on a geometric shape Geometry, you must first initialize two Projections, then create a CoordinateTransformation object and use it for transformation

sourceSR = osr.SpatialReference()

sourceSR.ImportFromEPSG(32612) #UTM 12N WGS84

targetSR = osr.SpatialReference()

targetSR.ImportFromEPSG(4326) #Geo WGS84

coordTrans = osr.CoordinateTransformation(sourceSR, targetSR)

geom.Transform(coordTrans)

But this code is very annoying ! It doesn't work in windows. There is a discussion in the foreigner's forum, saying that there is no problem in Linux, but Windows is dead, hehe. . .

There are a few more things to pay attention to:

Edit the Geometry at the appropriate time. It is best not to move it after the projection transformation.

To perform projection transformation on all Geometry in a DataSource, you have to do it one by one. Use a loop.

Writing your projection into a .prj file is actually very simple. First, MorphToESRI(), convert it into a string, then open a text file and write it in it. For example:

targetSR.MorphToESRI()

file = open('test.prj', 'w')

file.write(targetSR.ExportToWkt())

ffile.close()

The above is python gdal Tutorial: Geometry and projection content. For more related content, please pay attention to the PHP Chinese website (www.php.cn)!