python gdal tutorial: filters, simple spatial analysis, functions and modules

The Layer object has a method called SetAttributeFilter() that can filter out Features in the Layer that meet certain conditions. After setting the Filter, you can use the GetNextFeature() method to sequentially retrieve Features that meet the conditions. SetAttributeFilter(None) can clear a Filter. For example

>>> layer.GetFeatureCount()

42

>>> layer.SetAttributeFilter("cover = 'shrubs'")

>>> layer.GetFeatureCount() There are two types of spatial filters. One is SetSpatialFilter(), which filters a certain type of Feature. For example, filling in Polygon in the parameter is to select all Polygons in the Layer

There is also SetSpatialFilterRect(, , < ;maxx>, ), enter four coordinates as parameters, and you can select Feature

SetSpatialFilter(None) in the box to clear the spatial attribute filter.

For example, in the following code, layerAreas is polygon and layerSites is point

>>> featAreas = layerAreas.GetNextFeature()

>>> poly = featAreas.GetGeometryRef()

>> ; & gt; layersites.getfeatureCount ()

42

& gt; & gt; & gt; layersites.SetSpatialfilter (poly) & gt; & gt; & gt; & gt; Count () & gt; & gt; & gt; layersites.getfeatureCount ()

33

>>> layerSites.SetSpatialFilterRect(460000, 4590000, 490000, 4600000)

>>> layerSites.GetFeatureCount()

4

>>> layerSites.SetSpatialFilter( None)

>>> layerSites.GetFeatureCount()

42

There are also more complex Filters, such as executing SQL query statements ExecuteSQL(). With the powerful functions of SQL, more complex filters can be executed The task, such as the following code, is to select Features whose cover type is grass and sort them in descending order by ID number.

result = dsSites.ExecuteSQL("select * from sites where cover = 'grass' order by id desc")

resultFeat = result.GetNextFeature()

while resultFeat :

print resultFeat.GetField('id') print resultFeat.GetField('id')

       resultFeat = result.GetNextFeature()

  dsSites.ReleaseResultSet(result)

  42

  40

The last sentence ReleaseResultSet() is Release the query results and be sure to release them before executing the next SQL statement.

The following example counts the number of all Features whose cover is grass

>>> result = dsSites.ExecuteSQL("select count(*) from sites where cover = 'grass'")

> >> result.GetFeatureCount()

11

>>> result.GetFeature(0).GetField(0)

11

>>> dsSites.ReleaseResultSet(result)

List all different cover types

result = ds.ExecuteSQL("select distinct cover from sites") resultFeat = result.GetNextFeature()

while resultFeat:

print resultFeat.GetField(0)

resultFeat = result.GetNextFeature()

ds.ReleaseResultSet(result) shrubs

trees

rocks

grass

bare

Water

Count how many Features are there for each cover type

coverLayer = ds.ExecuteSQL ('select distinct cover from sites')

coverFeat = coverLayer.GetNextFeature()

while coverFeat:

cntLayer = ds.ExecuteSQL("select count(*) from sites where cover = ' " + coverFeat.GetField(0 ) + " ' ")

print coverFeat.GetField(0) + ' ' +print coverFeat.GetField(0) + ' ' + cntLayer.GetFeature(0).GetFieldAsString(0)

ds.ReleaseResultSet(cntLayer)

coverFeat = coverLayer.GetNextFeature()

ds.ReleaseResultSet(coverLayer)

shrubs 6

trees 11

rocks 6

grass 11

bare 6

water 2

Intersect determines whether two elements intersect

poly2.Intersect(poly1)

Returns 0 for disjoint and 1 for intersection.

Disjoint determines whether two elements are disjoint.
poly2.Disjoint(poly1)

Returns 1 for disjoint and 0 for intersection. Intersect is just the opposite

Touch means adjacent (edge ​​rubbing)

poly2.Touches(poly1)

Return 0 means no edge, return 1 means rub edge

Crosses cross, usually a line passes through a polygon

poly2.Crosses(line)

Returns 0 to indicate no crossing, returns 1 to indicate crossing

Within, one element is completely surrounded by another element

ptB.Within(poly1)

Returns 0 to indicate that the point is within Outside the polygon, return 1 to indicate that the point is within the polygon.

Contains contains, which is exactly the opposite of Within.

poly1.Contains(ptB)

Just change the main object and parameters.

Overlaps, it seems that there are only two polygons. The time required to overlap

poly2.Overlaps(poly3)

returns 0 to indicate non-overlap, and returns 1 to indicate overlap

Let’s look at simple geographical data processing geoprocessing

Polygon:

Intersection: poly3.Intersection(poly2)

Union: poly3.Union(poly2)

Difference: poly3.Difference(poly2)

Supplement: poly3.SymmetricDifference(poly2)

geometry:

.Buffer() to geometry Adding a buffer means turning the point line into a polygon, making it thicker

.Equal() Are the two geometries equal?

.Distance() Returns the shortest distance between two geometries

.GetEnvelope() Envelope, interesting, actually uses a box to enclose this geometric shape and returns four The coordinates of the corners (minx, maxx, miny, maxy)

Python’s function function, exception exception and module module can be found in any python textbook, so I won’t go into details here

The above is the python gdal tutorial: filtering , simple spatial analysis, functions and module content. For more related content, please pay attention to the PHP Chinese website (m.sbmmt.com)!