Commit 51de2c41 authored by Kris Vanneste's avatar Kris Vanneste
Browse files

Added new WCSConnection class.

Reimplemented __init__ method of WCSData in a more efficient way.
Added 'nodata_value', 'unit' and 'value_conversion' arguments to __init__ method of WCSData.
Added get_native_crs method to WCSData.
Renamed get_native_cellsize method of WCSData to get_native_cell_size, and simplified it.
Simplified get_native_bbox method of WCSData.
Replaced get_coverage method of WCSData with read_coverage.
Added get_url, close and __del__ methods to WCSData.
Renamed resx and resy property methods of WCSData to _resx and _resy.
Reipplemented to_gdal_raster_data method of WCSData.

git-svn-id: https://svn.seismo.oma.be/svn/seismo/mapping/layeredbasemap/trunk@7549 40b490c5-b4d9-47cb-8714-9bef99b524d5
parent 17e975ad
......@@ -26,7 +26,7 @@ from .polygon import PolygonData, MultiPolygonData
__all__ = ['GridData', 'MeshGridData', 'UnstructuredGridData',
'GdalRasterData', 'WCSData', 'MeshGridVectorData']
'GdalRasterData', 'WCSConnection', 'WCSData', 'MeshGridVectorData']
class GridData(BasemapData):
......@@ -1490,6 +1490,158 @@ class MeshGridVectorData(BasemapData):
return MeshGridVectorData(vx, vy, unit=unit)
class WCSConnection():
"""
WCS connection, providing information of available coverages (layers)
and of coverage metadata without actually having to read data from the
server
:param url:
str, base URL of WCS server
:param wcs_version:
str, WCS version
(default: '1.0.0')
:param verbose:
bool, whether or not to print information
(default: False)
"""
def __init__(self, url, wcs_version='1.0.0', verbose=False):
self.url = url
## Query server for some information
self.wcs_version = wcs_version
self.init_wcs(verbose=verbose)
def init_wcs(self, verbose=True):
"""
Initiate WCS connection
:param verbose:
bool, whether or not to print information
(default: True)
:return:
instance of :class:`owslib.wcs.WebCoverageService`
"""
from owslib.wcs import WebCoverageService
self.wcs = WebCoverageService(self.url, version=self.wcs_version)
if verbose:
print(self.get_layer_names())
def get_layer_names(self):
"""
Get list of available coverages (layers) on the server
:return:
list of strings
"""
return sorted(self.wcs.contents.keys())
def get_coverage(self, layer_name):
"""
Get a particular coverage (layer)
:param layer_name:
str, name of requested dataset on WCS server
return:
owslib coverage
"""
return self.wcs[layer_name]
get_layer = get_coverage
def get_coverage_info(self, layer_name):
"""
Get coordinate system and bounding box for a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
(crs, bbox, (dx, dy)) tuple
"""
crs = self.get_crs(layer_name)
bbox = self.get_bbox(layer_name)
(dx, dy) = self.get_cell_size(layer_name)
return (crs, bbox, (dx, dy))
def get_crs(self, layer_name):
"""
Get coordinate reference system of a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
instance of :class:`owslib.crs.Crs`
"""
coverage = self.get_coverage(layer_name)
crs = coverage.supportedCRS[0]
return crs
def get_srs(self, layer_name):
"""
Read Spatial Reference System of a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
osr SpatialReference object
"""
from mapping.geotools.coordtrans import get_epsg_srs
crs = self.get_crs(layer_name)
if crs.authority == 'EPSG':
srs = get_epsg_srs(crs.getcode())
return srs
def get_bbox(self, layer_name):
"""
Get bounding box (in native coordinates) of a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
(llx, lly, urx, ury) tuple in native coordinates
"""
coverage = self.get_coverage(layer_name)
bbox = coverage.boundingboxes[0]['bbox']
return bbox
def get_bbox_lonlat(self, layer_name):
"""
Get bounding box (in lon/lat) of a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
(llx, lly, urx, ury) tuple in degrees
"""
coverage = self.get_coverage(layer_name)
return coverage.boundingBoxWGS84
def get_cell_size(self, layer_name):
"""
Get native cell size of a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
:return:
(dx, dy) tuple
"""
coverage = self.get_coverage(layer_name)
dx = np.abs(float(coverage.grid.offsetvectors[0][0]))
dy = np.abs(float(coverage.grid.offsetvectors[1][1]))
return (dx, dy)
class WCSData(GdalRasterData):
"""
:param url:
......@@ -1514,100 +1666,68 @@ class WCSData(GdalRasterData):
list or tuple of floats: (lonmin, lonmax, latmin, latmax)
that will be used to determine bbox
(default: [])
:param nodata_value:
float, value indicating absence of data
(default: None)
:param unit:
str, measurement unit of gridded values
(default: "")
:param value_conversion:
function to apply to stored grid values.
Only applies if :param:`band_nr` is not null
(default: None)
:param wcs_version:
str, WCS version
(default: '1.0.0')
:param verbose:
bool, whether or not to print information
(default: True)
(default: False)
"""
def __init__(self, url, layer_name, down_sampling=1, resolution=None, band_nr=1, bbox=[], region=[],
wcs_version='1.0.0', verbose=True):
self.url = url
nodata_value=None, unit="", value_conversion=None,
wcs_version='1.0.0', verbose=False):
self.layer_name = layer_name
## Query server for some information
self.wcs_version = wcs_version
self.wcs = self.init_wcs(verbose=verbose)
self.wcs = WCSConnection(url, wcs_version=wcs_version, verbose=verbose)
if resolution:
dx, dy = self.get_native_cellsize()
dx, dy = self.get_native_cell_size()
if isinstance(resolution, (int, float)):
resolution = (resolution, resolution)
down_sampling = min(resolution[0] / dx, resolution[1] / dy)
self.down_sampling = max(1, down_sampling)
self.band_nr = band_nr
self._down_sampling = max(1, down_sampling)
## Set bounding box
if region:
bbox = self.get_bbox_from_region(region)
elif not bbox:
bbox = self.get_native_bbox()
self.bbox = self.add_margin_to_bbox(bbox)
self._bbox = self.add_margin_to_bbox(bbox)
## Read coverage
import urllib
response = self.get_coverage(self.layer_name)
if verbose:
if PY2:
unquote_func = urllib.unquote
else:
import urllib.parse
unquote_func = urllib.parse.unquote
print(unquote_func(response.geturl()))
#import tempfile
#fd = tempfile.NamedTemporaryFile(suffix=".tif")
#fd.write(response.read())
#fd.close()
#super(WCSData, self).__init__(fd.name, band_nr=band_nr, down_sampling=1)
mem_filename = self.read_coverage(verbose=verbose)
try:
super(WCSData, self).__init__(response.geturl(), band_nr=band_nr,
down_sampling=1)
## Old implementation using URL: very inefficient, as grid data
## will be read again from the server for each operation!
#super(WCSData, self).__init__(self._response.geturl(), band_nr=band_nr,
# down_sampling=1)
super(WCSData, self).__init__(mem_filename, band_nr=band_nr,
down_sampling=1, nodata_value=nodata_value,
unit=unit, value_conversion=value_conversion)
except:
print(response.read())
print(self._response.read())
raise
def init_wcs(self, verbose=True):
def get_native_crs(self):
"""
Initiate WCS connection
:param verbose:
bool, whether or not to print information
(default: True)
:return:
instance of :class:`owslib.wcs.WebCoverageService`
"""
from owslib.wcs import WebCoverageService
wcs = WebCoverageService(self.url, version=self.wcs_version)
if verbose:
print(sorted(wcs.contents.keys()))
return wcs
def get_coverage_info(self, layer_name):
"""
Get coordinate system and bounding box for a particular coverage
:param layer_name:
str, name of requested dataset on WCS server
Read native owslib coordinate reference system
:return:
(crs, bbox, (dx, dy)) tuple
instance of :class:`owslib.crs.Crs`
"""
coverage = self.wcs[layer_name]
crs = coverage.supportedCRS[0]
bbox = coverage.boundingboxes[0]['bbox']
dx = np.abs(float(coverage.grid.offsetvectors[0][0]))
dy = np.abs(float(coverage.grid.offsetvectors[1][1]))
return (crs, bbox, (dx, dy))
def get_crs(self):
coverage = self.wcs[self.layer_name]
crs = coverage.supportedCRS[0]
return crs
return self.wcs.get_crs(self.layer_name)
def get_native_srs(self):
"""
......@@ -1616,24 +1736,16 @@ class WCSData(GdalRasterData):
:return:
osr SpatialReference object
"""
from mapping.geotools.coordtrans import get_epsg_srs
return self.wcs.get_srs(self.layer_name)
crs = self.get_crs()
if crs.authority == 'EPSG':
srs = get_epsg_srs(crs.getcode())
return srs
def get_native_cellsize(self):
def get_native_cell_size(self):
"""
Read native cellsize
Read native cell size
:return:
(dx, dy) tuple
"""
coverage = self.wcs[self.layer_name]
dx = np.abs(float(coverage.grid.offsetvectors[0][0]))
dy = np.abs(float(coverage.grid.offsetvectors[1][1]))
return (dx, dy)
return self.wcs.get_cell_size(self.layer_name)
def get_native_bbox(self):
"""
......@@ -1642,9 +1754,7 @@ class WCSData(GdalRasterData):
:return:
(llx, lly, urx, ury) tuple in native coordinates
"""
coverage = self.wcs[self.layer_name]
bbox = coverage.boundingboxes[0]['bbox']
return bbox
return self.wcs.get_bbox(self.layer_name)
def get_bbox_from_region(self, region, margin_fraction=1./20):
"""
......@@ -1692,7 +1802,7 @@ class WCSData(GdalRasterData):
(llx, lly, urx, ury) tuple
"""
native_bbox = self.get_native_bbox()
resx, resy = self.resx, self.resy
resx, resy = self._resx, self._resy
llx, lly, urx, ury = bbox
llx = max(native_bbox[0], llx - resx)
urx = min(native_bbox[2], urx + resx)
......@@ -1700,65 +1810,90 @@ class WCSData(GdalRasterData):
ury = min(native_bbox[3], ury + resy)
return (llx, lly, urx, ury)
def get_coverage(self, layer_name):
def read_coverage(self, verbose=False):
"""
Fetch coverage from server
Fetch coverage from server and store as in-memory GeoTiff
:param layer_name:
str, name of requested dataset on WCS server
Note:
- relies on :prop:`_bbox`
- will store server response in :prop:`_response`
:return:
str, name of in-memory GeoTiff file
instance of :class:`owslib.util.ResponseWrapper`
"""
width, height = None, None
try:
crs = self.get_crs()
crs = self.get_native_crs()
except:
crs_code = "EPSG:31370"
else:
crs_code = crs.getcode()
format = "GeoTIFF"
return self.wcs.getCoverage(identifier=self.layer_name, width=width,
height=height, resx=self.resx, resy=self.resy, bbox=self.bbox,
format=format, crs=crs_code)
self._response = self.wcs.wcs.getCoverage(identifier=self.layer_name,
width=width, height=height,
resx=self._resx, resy=self._resy,
bbox=self._bbox,
format=format, crs=crs_code)
def to_gdal_raster_data(self, verbose=True):
"""
Convert to GDAL raster
if verbose:
import urllib
try:
## PY2
unquote_func = urllib.unquote
except:
## PY3
import urllib.parse
unquote_func = urllib.parse.unquote
print(unquote_func(self._response.geturl()))
:param verbose:
bool, whether or not to print some information
## Store fetched data in in-memory GeoTiff
mem_filename = "/vsimem/%s.tif" % self.layer_name
gdal.FileFromMemBuffer(mem_filename, self._response.read())
return mem_filename
def get_url(self):
"""
URL used to fetch coverage from server
:return:
instance of :class:`GdalRasterData`
str
"""
## Read coverage
import urllib
return self._response.geturl()
response = self.get_coverage(self.layer_name)
if verbose:
print(urllib.unquote(response.geturl()))
def to_gdal_raster_data(self, out_file):
"""
Convert to GDAL raster data and save to a GeoTiff file
#import tempfile
#fd = tempfile.NamedTemporaryFile(suffix=".tif")
#fd.write(response.read())
#fd.close()
#super(WCSData, self).__init__(fd.name, band_nr=band_nr, down_sampling=1)
:param out_file:
str, full path to output file (GeoTiff format)
try:
data = GdalRasterData(response.geturl(), band_nr=self.band_nr,
down_sampling=1)
except:
print(response.read())
raise
else:
return data
:return:
instance of :class:`GdalRasterData`
"""
with open(out_file, 'wb') as of:
of.write(self._response.read())
return GdalRasterData(out_file, band_nr=self.band_nr, down_sampling=1,
nodata_value=self.nodata_value, unit=self.unit,
value_conversion=self.value_conversion)
@property
def resx(self):
return self.get_native_cellsize()[0] * self.down_sampling
def _resx(self):
return self.get_native_cell_size()[0] * self._down_sampling
@property
def resy(self):
return self.get_native_cellsize()[1] * self.down_sampling
def _resy(self):
return self.get_native_cell_size()[1] * self._down_sampling
def close(self):
"""
Delete memory file
"""
gdal.Unlink(self.filespec)
def __del__(self):
## Delete memory file when instance is destroyed
self.close()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment