-
Notifications
You must be signed in to change notification settings - Fork 0
/
utils_wcs.py
188 lines (172 loc) · 6.7 KB
/
utils_wcs.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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
import numpy as np
from owslib.wcs import WebCoverageService
from owslib.util import Authentication
from shapely import wkt
# TODO Check how can it be improved
## TO READ WCS outputs
class WCS:
"""WCS object to get metadata etc and to get grid."""
def __init__(self, host, layer, username, password):
self.username = username
self.password = password
self.id = layer
self.wcs = (
WebCoverageService(
host,
version="1.0.0",
auth=Authentication(username=self.username, password=self.password),
)
if self.password and self.username
else WebCoverageService(host, version="1.0.0")
)
self.layer = self.wcs[self.id]
self.cx, self.cy = map(int, self.layer.grid.highlimits)
self.crs = self.layer.boundingboxes[0]["nativeSrs"]
self.bbox = self.layer.boundingboxes[0]["bbox"]
self.lx, self.ly, self.hx, self.hy = map(float, self.bbox)
self.resx, self.resy = (self.hx - self.lx) / self.cx, (
self.hy - self.ly
) / self.cy
self.width = self.cx
self.height = self.cy
def getw(self, fn):
"""Downloads raster and returns filename of written GEOTIFF in the tmp dir."""
gc = self.wcs.getCoverage(
identifier=self.id,
bbox=self.bbox,
format="GeoTIFF",
crs=self.crs,
width=self.width,
height=self.height,
)
f = open(fn, "wb")
f.write(gc.read())
f.close()
return fn
def getw_with_auth(self, fn):
"""Downloads raster and returns filename of written GEOTIFF in the tmp dir."""
gc = self.wcs.getCoverage(
identifier=self.id,
bbox=self.bbox,
format="GeoTIFF",
crs=self.crs,
width=self.width,
height=self.height,
auth=Authentication(username=self.username, password=self.password),
)
f = open(fn, "wb")
f.write(gc.read())
f.close()
return fn
## TO handle transects
class LS:
"""Intersection on grid line"""
def __init__(self, awkt, crs, host, layer, username, password, sampling=1):
self.wwkt = awkt
self.crs = crs
self.gs = WCS(
host, layer, username, password
) # Initiates WCS service to get some parameters about the grid.
self.sampling = sampling
def line(self):
"""Creates WCS parameters and sample coordinates for cells in raster based on line input."""
self.ls = wkt.loads(self.wwkt)
self.ax, self.ay, self.bx, self.by = self.ls.bounds
# TODO http://stackoverflow.com/questions/13439357/extract-point-from-raster-in-gdal
"""if first x is larger than second, coordinates will be flipped during process of defining bounding box !!!!
next lines introduce a boolean flip variable used in the last part of this process"""
flipx = False
flipy = False
ax, bx = self.ls.coords.xy[0]
ay, by = self.ls.coords.xy[1]
if ax >= bx:
flipx = True
if ay >= by:
flipy = True
"""get raster coordinates"""
self.ax = (
self.ax - self.gs.lx
) # coordinates minus coordinates of raster, start from 0,0
self.ay = self.ay - self.gs.ly
self.bx = self.bx - self.gs.lx
self.by = self.by - self.gs.ly
self.x1, self.y1 = int(self.ax // self.gs.resx), int(self.ay // self.gs.resy)
self.x2, self.y2 = (
int(self.bx // self.gs.resx) + 1,
int(self.by // self.gs.resy) + 1,
)
self.gs.bbox = (
self.x1 * self.gs.resx + self.gs.lx,
self.y1 * self.gs.resy + self.gs.ly,
self.x2 * self.gs.resx + self.gs.lx,
self.y2 * self.gs.resy + self.gs.ly,
)
self.gs.width = abs(self.x2 - self.x1) # difference of x cells
self.gs.height = abs(self.y2 - self.y1)
""" here we go back to our line again instead of calculating bbox for wcs request."""
self.ax, self.bx = self.ls.coords.xy[0]
self.ay, self.by = self.ls.coords.xy[1]
# coordinates minus coordinates of raster, start from 0,0
self.ax = self.ax - self.gs.lx
self.ay = self.ay - self.gs.ly
self.bx = self.bx - self.gs.lx
self.by = self.by - self.gs.ly
if flipx and flipy: # who draws these lines?
# top right to bottom left
self.x2, self.y2 = int(self.bx // self.gs.resx), int(
self.by // self.gs.resy
)
self.x1, self.y1 = (
int(self.ax // self.gs.resx) + 1,
int(self.ay // self.gs.resy) + 1,
)
elif flipx:
# bottom right to top left
self.x2, self.y1 = int(self.bx // self.gs.resx), int(
self.ay // self.gs.resy
)
self.x1, self.y2 = (
int(self.ax // self.gs.resx) + 1,
int(self.by // self.gs.resy) + 1,
)
elif flipy:
# top left to bottom right
self.x1, self.y2 = int(self.ax // self.gs.resx), int(
self.by // self.gs.resy
)
self.x2, self.y1 = (
int(self.bx // self.gs.resx) + 1,
int(self.ay // self.gs.resy) + 1,
)
else:
# normal
self.x1, self.y1 = int(self.ax // self.gs.resx), int(
self.ay // self.gs.resy
)
self.x2, self.y2 = (
int(self.bx // self.gs.resx) + 1,
int(self.by // self.gs.resy) + 1,
)
# From upperright to lower left x values become negative
# Subdivide the line into sampling points of the raster.
# Takes longest dimension and uses number of cells * sampling as the
# number of subdivisions.
# Grid of subdivions is pixel grid - 0.5
self.subdiv = int(max(self.gs.width, self.gs.height)) * self.sampling
self.xlist = np.linspace(
(self.ax / self.gs.resx) - min(self.x1, self.x2),
(self.bx / self.gs.resx) - min(self.x1, self.x2),
num=self.subdiv,
)
self.ylist = np.linspace(
(self.ay / self.gs.resy) - min(self.y1, self.y2),
(self.by / self.gs.resy) - min(self.y1, self.y2),
num=self.subdiv,
)
def getraster(self, fname, all_box=False):
"""Returns values of line intersection on downlaoded geotiff from wcs."""
if self.gs.username and self.gs.password:
self.gs.getw_with_auth(fname)
else:
self.gs.getw(fname)
return