### Sidebar

spatial-ecology.org

Trainings:

Learn:

Data:

Community:
teachers
students
projects
Matera 2015
Vancouver 2015
Santa Barbara 2015
Site design
Hands on:
Installations

Donations
USD

GBP

wiki:python:basicpython5

# Using GDAL/OGR API for GIS analysis

Material prepared by Stephen Weston

```cd ~/ost4sem/exercise/basic_adv_gdalogr/fagus_sylvatica
wget  https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/python/gdal_example.py
wget  https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/python/ogr_example.py
wget  https://raw.githubusercontent.com/selvaje/spatial-ecology-codes/master/wiki/python/polygon_means.py```

### Reading raster images and return Min Max Mean

wiki/python/gdal_example.py
```#!/us/bin/env python
import os, sys
from osgeo import gdal
from osgeo.gdalconst import *

driver = gdal.GetDriverByName('GTiff')
driver.Register()

'2020_TSSP_IP-ENS-A2-SP20_43023435.tif')
if ds is None:
print 'Error opening', fn
sys.exit(1)

cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
print cols, rows, bands

geotransform = ds.GetGeoTransform()
print geotransform

for i in range(bands):
band = ds.GetRasterBand(i + 1)
d = band.ReadAsArray(0, 0, band.XSize, band.YSize)
print d.shape, d.dtype
print d.min(), d.max(), d.mean()```

### Reading a vector file and return a coordinates of the centroid

wiki/python/ogr_example.py
```#!/usr/bin/env python
import os, sys
from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(fn, 0)
if dataSource is None:
print 'Error opening', fn
sys.exit(1)

layer = dataSource.GetLayer(0)
numFeatures = layer.GetFeatureCount()
print 'Feature count:', numFeatures

extend = layer.GetExtent()
print 'Extend:', extend

feature = layer.GetNextFeature()
while feature:
print feature.GetField('id')
geometry = feature.GetGeometryRef()
print geometry.Centroid()
print geometry.ExportToWkt()
feature.Destroy()
feature = layer.GetNextFeature()

dataSource.Destroy()```

### Demonstrate how to create a NumPy mask array

```#!/usr/bin/env python
import numpy as np
import numpy.ma as ma

x = np.array([[1, 2, 3], [3, -1, 5]])
print x.mean()

mask = np.array([[0, 0, 0], [0, 1, 0]])
print mx.mean()
print mx

y = np.array([1, 2, 3, 3, 5])
print y.mean()

print x[x >= 0].mean()

####
mx = ma.masked_where(x < 0, x)
print mx.mean()

mmx = ma.masked_where(mx > 2, mx)
print mmx

print mmx.fill_value```

### Compute statistic inside polygons

`gdal_rasterize -te 2632000.000 1420000.000 5972000.000 5450000.00 -tap -clump -tr 1000 1000 -l poly_fr poly_fr.shp poly_fr.tif`
wiki/python/polygon_means.py
```#!/us/bin/env python
#
# Rasterize shape file using gdal_rasterize:
#
# \$ gdal_rasterize -te 2632000.000 1420000.000 5972000.000 5450000.000 \
#    -tap -clump -tr 1000 1000 -l poly_fr poly_fr.shp poly_fr.tif
#
import os, sys
from osgeo import gdal
from osgeo.gdalconst import *
import numpy as np
import numpy.ma as ma

driver = gdal.GetDriverByName('GTiff')
driver.Register()

if polyds is None:
print 'Error opening', rasterized
sys.exit(1)

'2020_TSSP_IP-ENS-A2-SP20_43023435.tif')
if ds is None:
print 'Error opening', raster
sys.exit(1)

# Get raster data
band = ds.GetRasterBand(1)
rasterD = band.ReadAsArray(0, 0, band.XSize, band.YSize)

# Get shape data
band = polyds.GetRasterBand(1)
shapes = band.ReadAsArray(0, 0, band.XSize, band.YSize)

print "total value count:", band.XSize * band.YSize
print "ignored value count:", (shapes == 0.0).sum()

# Iterate over the different polygons
for d in np.unique(shapes):
print 'Polygon:', d
# if d == 0.0: continue