-
Notifications
You must be signed in to change notification settings - Fork 0
/
zonal-stats-benchmark.py
60 lines (44 loc) · 1.55 KB
/
zonal-stats-benchmark.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# -*- coding: utf-8 -*-
__author__ = 'Douglas Uba'
__email__ = 'douglas.uba@inpe.br'
import matplotlib.pyplot as plt
from utils import create_image, create_random_circles, Timer
import numpy as np
from affine import Affine
from rasterstats import zonal_stats
# User-defined statistics example
def model(x):
value = np.max(x) - np.min(x)
if np.ma.is_masked(value):
value = None
return value
def compute_stats(image, geoms):
# Get Affine object in order to run zonal_stats
aff = Affine.from_gdal(*image.GetGeoTransform())
# Extract values
values = image.ReadAsArray()
# Get no-data value
nodata = image.GetRasterBand(1).GetNoDataValue()
# Create WKT representation for each polygon
wkts = []
for g in geoms:
wkts.append(g.ExportToWkt())
# Compute stats for each polygon
stats = zonal_stats(wkts, values, stats=['min', 'max', 'count', 'std'], affine=aff,
nodata=nodata, raster_out=False, prefix='', add_stats={'model' : model}, geojson_out=True)
return stats
# Define geo extent
llx = -75.0; lly = -35.50; urx = -34.0; ury = 5.54
extent = [llx, lly, urx, ury]
# Define resolution
resolution = 2.0 # km
# Create image
image = create_image(extent, resolution, driver='MEM')
# Create random geometries
geoms = create_random_circles(extent, min=0.5, max=30.0, n=32000)
with Timer():
stats = compute_stats(image, geoms)
import json
result = {'type': 'FeatureCollection', 'features': stats}
with open('./data/geoms-benchmark.geojson', 'w') as file:
json.dump(result, file)