77
88from uuid import UUID
99from logging import debug
10- from io import StringIO , BytesIO
10+ from io import BytesIO
1111import urllib .parse
1212import json
1313
1414import requests as req
1515import geopandas as gpd
16- import matplotlib .pyplot as plt
17- import cartopy .crs as ccrs
1816from owslib .util import Authentication
1917from owslib .wcs import WebCoverageService
20- from owslib .wms import WebMapService
21- import rasterio
18+ import rasterio .io
2219from vega import VegaLite
2320import numpy as np
2421from PIL import Image
22+ import xarray as xr
2523
2624from geoengine .types import InternalDatasetId , ProvenanceOutput , QueryRectangle , ResultDescriptor
2725from geoengine .auth import get_session
2826from geoengine .error import GeoEngineException , MethodNotCalledOnPlotException , MethodNotCalledOnRasterException , \
29- MethodNotCalledOnVectorException , SpatialReferenceMismatchException
27+ MethodNotCalledOnVectorException , SpatialReferenceMismatchException , check_response_for_error
3028from geoengine .datasets import StoredDataset , UploadId
3129
3230
@@ -120,15 +118,15 @@ def __get_wfs_url(self, bbox: QueryRectangle) -> str:
120118 version = "2.0.0" ,
121119 request = 'GetFeature' ,
122120 outputFormat = 'application/json' ,
123- typeNames = f'registry: { self .__workflow_id } ' ,
121+ typeNames = f'{ self .__workflow_id } ' ,
124122 bbox = bbox .bbox_str ,
125123 time = bbox .time_str ,
126124 srsName = bbox .srs ,
127125 queryResolution = f'{ bbox .resolution [0 ]} ,{ bbox .resolution [1 ]} '
128126 )
129127
130128 wfs_url = req .Request (
131- 'GET' , url = f'{ session .server_url } /wfs' , params = params ).prepare ().url
129+ 'GET' , url = f'{ session .server_url } /wfs/ { self . __workflow_id } ' , params = params ).prepare ().url
132130
133131 debug (f'WFS URL:\n { wfs_url } ' )
134132
@@ -165,17 +163,20 @@ def get_dataframe(self, bbox: QueryRectangle) -> gpd.GeoDataFrame:
165163
166164 data_response = req .get (wfs_url , headers = session .auth_header )
167165
168- def geo_json_with_time_to_geopandas (data_response ):
166+ check_response_for_error (data_response )
167+
168+ data = data_response .json ()
169+
170+ def geo_json_with_time_to_geopandas (geo_json ):
169171 '''
170172 GeoJson has no standard for time, so we parse the when field
171173 separately and attach it to the data frame as columns `start`
172174 and `end`.
173175 '''
174176
175- data = gpd .read_file ( StringIO ( data_response . text ) )
177+ data = gpd .GeoDataFrame . from_features ( geo_json )
176178 data = data .set_crs (bbox .srs , allow_override = True )
177179
178- geo_json = data_response .json ()
179180 start = [f ['when' ]['start' ] for f in geo_json ['features' ]]
180181 end = [f ['when' ]['end' ] for f in geo_json ['features' ]]
181182
@@ -186,65 +187,16 @@ def geo_json_with_time_to_geopandas(data_response):
186187
187188 return data
188189
189- return geo_json_with_time_to_geopandas (data_response )
190-
191- def plot_image (self , bbox : QueryRectangle , ax : plt .Axes = None , timeout = 3600 ) -> plt .Axes :
192- '''
193- Query a workflow and return the WMS result as a matplotlib image
194-
195- Params:
196- timeout - - HTTP request timeout in seconds
197- '''
198-
199- if not self .__result_descriptor .is_raster_result ():
200- raise MethodNotCalledOnRasterException ()
201-
202- session = get_session ()
203-
204- wms_url = f'{ session .server_url } /wms'
205-
206- def srs_to_projection (srs : str ) -> ccrs .Projection :
207- fallback = ccrs .PlateCarree ()
208-
209- [authority , code ] = srs .split (':' )
210-
211- if authority != 'EPSG' :
212- return fallback
213- try :
214- return ccrs .epsg (code )
215- except ValueError :
216- return fallback
217-
218- if ax is None :
219- ax = plt .axes (projection = srs_to_projection (bbox .srs ))
220-
221- wms = WebMapService (wms_url ,
222- version = '1.3.0' ,
223- xml = self .__faux_capabilities (wms_url , str (self ), bbox ),
224- auth = Authentication (auth_delegate = session .requests_bearer_auth ()),
225- timeout = timeout )
226-
227- # TODO: incorporate spatial resolution (?)
228-
229- ax .add_wms (wms ,
230- layers = [str (self )],
231- wms_kwargs = {
232- 'time' : urllib .parse .quote (bbox .time_str ),
233- # 'bbox': bbox.bbox_str
234- 'crs' : bbox .srs
235- })
236-
237- ax .set_xlim (bbox .xmin , bbox .xmax )
238- ax .set_ylim (bbox .ymin , bbox .ymax )
239-
240- return ax
190+ return geo_json_with_time_to_geopandas (data )
241191
242192 def wms_get_map_as_image (self , bbox : QueryRectangle , colorizer_min_max : Tuple [float , float ] = None ) -> Image :
243193 '''Return the result of a WMS request as a PIL Image'''
244194
245195 wms_request = self .__wms_get_map_request (bbox , colorizer_min_max )
246196 response = req .Session ().send (wms_request )
247197
198+ check_response_for_error (response )
199+
248200 return Image .open (BytesIO (response .content ))
249201
250202 def __wms_get_map_request (self ,
@@ -291,7 +243,7 @@ def __wms_get_map_request(self,
291243
292244 return req .Request (
293245 'GET' ,
294- url = f'{ session .server_url } /wms' ,
246+ url = f'{ session .server_url } /wms/ { str ( self ) } ' ,
295247 params = params ,
296248 headers = session .auth_header
297249 ).prepare ()
@@ -306,56 +258,6 @@ def wms_get_map_curl(self, bbox: QueryRectangle, colorizer_min_max: Tuple[float,
306258 headers = " -H " .join (headers )
307259 return command .format (method = wms_request .method , headers = headers , uri = wms_request .url )
308260
309- @classmethod
310- def __faux_capabilities (cls , wms_url : str , layer_name : str , bbox : QueryRectangle ) -> str :
311- '''Create an XML file with faux capabilities to list the layer with `layer_name`'''
312-
313- return '''
314- <WMS_Capabilities xmlns="http://www.opengis.net/wms" xmlns:sld="http://www.opengis.net/sld" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/wms http://schemas.opengis.net/wms/1.3.0/capabilities_1_3_0.xsd http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/sld_capabilities.xsd">
315- <Service>
316- <Name>WMS</Name>
317- <Title>Geo Engine WMS</Title>
318- </Service>
319- <Capability>
320- <Request>
321- <GetCapabilities>
322- <Format>text/xml</Format>
323- <DCPType>
324- <HTTP>
325- <Get>
326- <OnlineResource xlink:href="{wms_url}"/>
327- </Get>
328- </HTTP>
329- </DCPType>
330- </GetCapabilities>
331- <GetMap>
332- <Format>image/png</Format>
333- <DCPType>
334- <HTTP>
335- <Get>
336- <OnlineResource xlink:href="{wms_url}"/>
337- </Get>
338- </HTTP>
339- </DCPType>
340- </GetMap>
341- </Request>
342- <Exception>
343- <Format>XML</Format>
344- <Format>INIMAGE</Format>
345- <Format>BLANK</Format>
346- </Exception>
347- <Layer queryable="1">
348- <Name>{layer_name}</Name>
349- <Title>{layer_name}</Title>
350- <CRS>{crs}</CRS>
351- <BoundingBox CRS="EPSG:4326" minx="-90.0" miny="-180.0" maxx="90.0" maxy="180.0"/>
352- </Layer>
353- </Capability>
354- </WMS_Capabilities>
355- ''' .format (wms_url = wms_url ,
356- layer_name = layer_name ,
357- crs = bbox .srs )
358-
359261 def plot_chart (self , bbox : QueryRectangle ) -> VegaLite :
360262 '''
361263 Query a workflow and return the plot chart result as a vega plot
@@ -372,18 +274,24 @@ def plot_chart(self, bbox: QueryRectangle) -> VegaLite:
372274
373275 plot_url = f'{ session .server_url } /plot/{ self } ?bbox={ spatial_bounds } &time={ time } &spatialResolution={ resolution } '
374276
375- response = req .get (plot_url , headers = session .auth_header ).json ()
277+ response = req .get (plot_url , headers = session .auth_header )
278+
279+ check_response_for_error (response )
280+
281+ response = response .json ()
376282
377283 vega_spec = json .loads (response ['data' ]['vegaString' ])
378284
379285 return VegaLite (vega_spec )
380286
381- def get_array (self , bbox : QueryRectangle , timeout = 3600 ) -> np . ndarray :
287+ def __get_wcs_tiff_as_memory_file (self , bbox : QueryRectangle , timeout = 3600 ) -> rasterio . io . MemoryFile :
382288 '''
383- Query a workflow and return the raster result as a numpy array
289+ Query a workflow and return the raster result as a memory mapped GeoTiff
384290
385- Params:
386- timeout - - HTTP request timeout in seconds
291+ Parameters
292+ ----------
293+ bbox : A bounding box for the query
294+ timeout : HTTP request timeout in seconds
387295 '''
388296
389297 if not self .__result_descriptor .is_raster_result ():
@@ -412,11 +320,47 @@ def get_array(self, bbox: QueryRectangle, timeout=3600) -> np.ndarray:
412320 resx = resx ,
413321 resy = resy ,
414322 timeout = timeout ,
415- )
323+ ).read ()
324+
325+ # response is checked via `raise_on_error` in `getCoverage` / `openUrl`
326+
327+ memory_file = rasterio .io .MemoryFile (response )
328+
329+ return memory_file
330+
331+ def get_array (self , bbox : QueryRectangle , timeout = 3600 ) -> np .ndarray :
332+ '''
333+ Query a workflow and return the raster result as a numpy array
334+
335+ Parameters
336+ ----------
337+ bbox : A bounding box for the query
338+ timeout : HTTP request timeout in seconds
339+ '''
340+
341+ with self .__get_wcs_tiff_as_memory_file (bbox , timeout ) as memfile , memfile .open () as dataset :
342+ array = dataset .read (1 )
416343
417- with rasterio .io .MemoryFile (response .read ()) as memfile , memfile .open () as dataset :
418344 # TODO: map nodata values to NaN?
419- return dataset .read (1 )
345+
346+ return array
347+
348+ def get_xarray (self , bbox : QueryRectangle , timeout = 3600 ) -> np .ndarray :
349+ '''
350+ Query a workflow and return the raster result as a georeferenced xarray
351+
352+ Parameters
353+ ----------
354+ bbox : A bounding box for the query
355+ timeout : HTTP request timeout in seconds
356+ '''
357+
358+ with self .__get_wcs_tiff_as_memory_file (bbox , timeout ) as memfile , memfile .open () as dataset :
359+ data_array = xr .open_rasterio (dataset )
360+
361+ # TODO: add time information to dataset
362+
363+ return data_array .load ()
420364
421365 def get_provenance (self ) -> List [ProvenanceOutput ]:
422366 '''
@@ -457,10 +401,11 @@ def save_as_dataset(self, bbox: QueryRectangle, name: str, description: str = ''
457401 url = f'{ session .server_url } /datasetFromWorkflow/{ self .__workflow_id } ' ,
458402 json = request_body ,
459403 headers = session .auth_header ,
460- ). json ()
404+ )
461405
462- if 'error' in response :
463- raise GeoEngineException (response )
406+ check_response_for_error (response )
407+
408+ response = response .json ()
464409
465410 return StoredDataset (
466411 dataset_id = InternalDatasetId .from_response (response ['dataset' ]),
0 commit comments