-
Notifications
You must be signed in to change notification settings - Fork 8
/
analysis.py
78 lines (61 loc) · 2.2 KB
/
analysis.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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import numpy as np
import rasterio
import matplotlib.pyplot as plt
from itertools import product
class Quadrats(object):
def __init__(self, filename=None, data=None, band=1, dx=None, dy=None):
if filename is not None:
self.load_data(filename, band=band)
elif data is not None:
self.data = data
else:
self.data = None
if dx is not None:
self.make_quadrats(dx, dy)
else:
self.quadrats = []
def load_data(self, filename, band=1):
with rasterio.open(filename, 'r') as src:
self.data = src.read(band)
def make_quadrats(self, dx, dy=None):
if dy is None:
dy = dx
ny, nx = self.data.shape
rows = np.arange(0, ny-dy+1, step=dy)
cols = np.arange(0, nx-dx+1, step=dx)
# TODO: Implement offset quadrats
# rows += row_offset
# cols += col_offset
# rows = rows[rows <= ny]
# cols = cols[cols <= nx]
self.quadrats = [self.data[r:r+dy, c:c+dx] for r in rows for c in cols]
def map_quadrats(self, func, **kwargs):
return [func(q, **kwargs) for q in self.quadrats]
def plot(self, values, ax=None, **kwargs):
ny, nx = self.data.shape
dy, dx = self.quadrats[0].shape
rows = np.arange(0, ny-dy+1, step=dy) + dy / 2
cols = np.arange(0, nx-dx+1, step=dx) + dx / 2
y = [p[0] for p in product(rows, cols)]
x = [p[1] for p in product(rows, cols)]
if ax is None:
plt.figure()
ax = plt.gca()
ax.scatter(x, y, c=values, **kwargs)
ax.set_xlim([0, nx])
ax.set_ylim([0, ny])
ax.invert_yaxis()
def quiver(self, u, v, ax=None, **kwargs):
ny, nx = self.data.shape
dy, dx = self.quadrats[0].shape
rows = np.arange(0, ny-dy+1, step=dy) + dy / 2
cols = np.arange(0, nx-dx+1, step=dx) + dx / 2
y = [p[0] for p in product(rows, cols)]
x = [p[1] for p in product(rows, cols)]
if ax is None:
plt.figure()
ax = plt.gca()
ax.quiver(x, y, u, v, **kwargs)
ax.set_xlim([0, nx])
ax.set_ylim([0, ny])
ax.invert_yaxis()