forked from mdbartos/RIPS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rect_grid.py
116 lines (93 loc) · 4.2 KB
/
rect_grid.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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import numpy as np
import pandas as pd
from shapely import geometry
import geopandas as gpd
def rect_grid(bbox, hdim, vdim=None, out='polygon', pointref='centroid',
how='fixed_step', anchor_point='ll', endpoint=True):
if isinstance(bbox, (gpd.geoseries.GeoSeries,
gpd.geodataframe.GeoDataFrame)):
bbox = bbox.total_bounds
elif isinstance(bbox, (geometry.Point, geometry.LineString,
geometry.LinearRing, geometry.Polygon,
geometry.MultiPoint, geometry.MultiLineString,
geometry.MultiPolygon)):
bbox = bbox.bounds
x0 = bbox[0]
y0 = bbox[1]
x1 = bbox[2]
y1 = bbox[3]
if how == 'fixed_step':
hstep = hdim
if vdim == None:
vstep = hstep
else:
vstep = vdim
if anchor_point.lower()[1] == 'l':
x = np.arange(x0, x1, hstep, dtype=float)
elif anchor_point.lower()[1] == 'r':
x = np.arange(x1, x0, -hstep, dtype=float)[::-1]
if anchor_point.lower()[0] == 'l':
y = np.arange(y0, y1, vstep, dtype=float)
elif anchor_point.lower()[0] == 'u':
y = np.arange(y1, y0, -vstep, dtype=float)[::-1]
elif how == 'fixed_number':
if vdim == None:
vdim = hdim
if anchor_point.lower()[1] == 'l':
x, hstep = np.linspace(x0, x1, hdim, retstep=True, dtype=float, endpoint=endpoint)
elif anchor_point.lower()[1] == 'r':
x, hstep = np.linspace(x1, x0, hdim, retstep=True, dtype=float, endpoint=endpoint)
x, hstep = x[::-1], -hstep
if anchor_point.lower()[0] == 'l':
y, vstep = np.linspace(y0, y1, vdim, retstep=True, dtype=float, endpoint=endpoint)
elif anchor_point.lower()[0] == 'u':
y, vstep = np.linspace(y1, y0, vdim, retstep=True, dtype=float, endpoint=endpoint)
y, vstep = y[::-1], -vstep
xy_ll = np.vstack(np.dstack(np.meshgrid(x, y)))
if out == 'point':
if pointref == 'centroid':
out_arr = np.column_stack([xy_ll[:,0] + hstep/2.0,
xy_ll[:,1] + vstep/2.0])
return gpd.GeoSeries(map(geometry.asPoint, out_arr))
elif pointref == 'll':
return gpd.GeoSeries(map(geometry.asPoint, xy_ll))
elif pointref == 'lr':
out_arr = np.column_stack([xy_ll[:,0] + hstep, xy_ll[:,1]])
return gpd.GeoSeries(map(geometry.asPoint, out_arr))
elif pointref == 'ur':
out_arr = np.column_stack([xy_ll[:,0] + hstep, xy_ll[:,1] + vstep])
return gpd.GeoSeries(map(geometry.asPoint, out_arr))
elif pointref == 'ul':
out_arr = np.column_stack([xy_ll[:,0], xy_ll[:,1] + vstep])
return gpd.GeoSeries(map(geometry.asPoint, out_arr))
elif out == 'line':
# if how == 'fixed_step':
# y1 = y1 + vstep
# x1 = x1 + hstep
vlines = np.hstack([
np.column_stack([x, np.repeat(y0, len(x))]),
np.column_stack([x, np.repeat(y1, len(x))])
])
vlines = vlines.reshape(vlines.shape[0], 2, 2)
hlines = np.hstack([
np.column_stack([np.repeat(x0, len(y)), y]),
np.column_stack([np.repeat(x1, len(y)), y])
])
hlines = hlines.reshape(hlines.shape[0], 2, 2)
out_arr = np.vstack([vlines, hlines])
del vlines
del hlines
return gpd.GeoSeries(map(geometry.asLineString, out_arr))
elif out == 'polygon':
xy_ur = np.column_stack([xy_ll[:,0] + hstep, xy_ll[:,1] + vstep])
out_arr = np.column_stack([
xy_ur,
np.column_stack([xy_ll[:,0], xy_ll[:,1] + vstep]),
xy_ll,
np.column_stack([xy_ll[:,0] + hstep, xy_ll[:,1]]),
xy_ur
])
del xy_ll
del xy_ur
out_arr = out_arr.reshape(out_arr.shape[0], out_arr.shape[1]/2, 2)
return gpd.GeoSeries(map(geometry.asPolygon, out_arr))