User Tools

Site Tools


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/masked.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()
 
fn = os.path.join(os.environ.get('DATADIR', '.'),
                  '2020_TSSP_IP-ENS-A2-SP20_43023435.tif')
ds = gdal.Open(fn, GA_ReadOnly)
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')
fn = os.path.join(os.environ.get('DATADIR', '.'), 'poly_fr.shp')
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

wiki/python/masked.py
#!/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]])
mx = ma.masked_array(x, mask=mask)
print mx.mean()
print mx
 
y = np.array([1, 2, 3, 3, 5])
print y.mean()
 
print x[x >= 0].mean()
print x[mask == 0].mean()
 
####
mx = ma.masked_where(x < 0, x)
print mx.mean()
 
# Mask a masked array
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()
 
rasterized = os.path.join(os.environ.get('DATADIR', '.'), 'poly_fr.tif')
polyds = gdal.Open(rasterized, GA_ReadOnly)
if polyds is None:
    print 'Error opening', rasterized
    sys.exit(1)
 
raster = os.path.join(os.environ.get('DATADIR', '.'),
                      '2020_TSSP_IP-ENS-A2-SP20_43023435.tif')
ds = gdal.Open(raster, GA_ReadOnly)
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
    mask = shapes != d
    mx = ma.masked_array(rasterD, mask=mask)
    print (shapes == d).sum()
    print mx.count()
    print mx.mean()
wiki/python/basicpython5.txt · Last modified: 2020/07/17 06:39 (external edit)