Commit f2cf6ec5 authored by Kris Vanneste's avatar Kris Vanneste
Browse files

Added missing docstrings to MeshGridData and GdalRasterData methods.

Implemented reproject_coordinates method of MeshGridData.
Updated extract_contour_lines method of MeshGridData to more recent versions of matplotlib.
Use center_lons/lats instead of lons/lats in to_multi_point_data and extract_contour_intervals methods of MeshGridData.
Improved export method of GdalRasterData for grids with bounding box and downsampling applied.

git-svn-id: https://svn.seismo.oma.be/svn/seismo/mapping/layeredbasemap/trunk@7462 40b490c5-b4d9-47cb-8714-9bef99b524d5
parent baa38cea
......@@ -341,18 +341,36 @@ class MeshGridData(GridData):
if apply_mask and self.is_masked():
mask = self.values.mask
values = self.values[~mask].data
lons = self.lons[~mask]
lats = self.lats[~mask]
lons = self.center_lons[~mask]
lats = self.center_lats[~mask]
else:
values = self.values.flatten()
lons = self.lons.flatten()
lats = self.lats.flatten()
lons = self.center_lons.flatten()
lats = self.center_lats.flatten()
return MultiPointData(lons, lats, list(values))
def mask_oceans(self, resolution, mask_lakes=False, grid_spacing=1.25):
"""
Apply ocean mask to grid
:param resolution:
char, basemap resolution ('c', 'l', 'i', 'h' or 'f')
:param mask_lakes:
bool, whether or not to mask inland water bodies too
(default: False)
:param grid_spacing:
float, land/see mask grid spacing in minutes (10, 5, 2.5 or 1.25)
(default: 1.25)
:return:
instance of :class:`MeshGridData`
"""
from mpl_toolkits.basemap import maskoceans
masked_values = maskoceans(self.lons, self.lats, self.values, inlands=mask_lakes, resolution=resolution, grid=grid_spacing)
return MeshGridData(self.lons, self.lats, masked_values)
masked_values = maskoceans(self.center_lons, self.center_lats,
self.values, inlands=mask_lakes,
resolution=resolution, grid=grid_spacing)
return MeshGridData(self.center_lons, self.center_lats, masked_values,
unit=self.unit)
def is_masked(self):
"""
......@@ -445,8 +463,23 @@ class MeshGridData(GridData):
self.apply_mask(mask)
def reproject(self, target_srs):
pass
def reproject_coordinates(self, target_srs):
"""
Reproject grid coordinates
:param target_srs:
instance of :class:`osgeo.osr.SpatialReference`
target spatial reference system
:return:
(XX, YY) tuple of 2D arrays, reprojected grid coordinates
"""
from mapping.geotools.coordtrans import transform_mesh_coordinates
X, Y = transform_mesh_coordinates(self.srs, target_srs,
self.center_lons, self.center_lats)
return X, Y
def interpolate(self, xout, yout, srs=None, checkbounds=False, masked=True,
order=1):
......@@ -560,19 +593,38 @@ class MeshGridData(GridData):
:return:
list with instances of :class:`MultiLineData`
"""
import matplotlib._cntr as cntr
import matplotlib
contour_engine = cntr.Cntr(self.lons, self.lats, self.values)
contour_lines = []
for level in levels:
nlist = contour_engine.trace(level, level, 0)
nseg = len(nlist) // 2
segs = nlist[:nseg]
contour_line = MultiLineData([], [])
for seg in segs:
cl = LineData(seg[:,0], seg[:,1], value=level)
contour_line.append(cl)
contour_lines.append(contour_line)
if matplotlib.__version__ <= '2.1.0':
import matplotlib._cntr as cntr
contour_engine = cntr.Cntr(self.center_lons, self.center_lats, self.values)
for level in levels:
nlist = contour_engine.trace(level, level, 0)
nseg = len(nlist) // 2
segs = nlist[:nseg]
contour_line = MultiLineData([], [])
for seg in segs:
cl = LineData(seg[:,0], seg[:,1], value=level)
contour_line.append(cl)
contour_lines.append(contour_line)
else:
import numpy.ma as ma
from matplotlib._contour import QuadContourGenerator
values = ma.asarray(self.values)
contour_engine = QuadContourGenerator(self.center_lons, self.center_lats,
values.filled(), None, False, 0)
for level in levels:
segments = contour_engine.create_contour(level)
contour_line = MultiLineData([], [])
for seg in segments:
cl = LineData(seg[:,0], seg[:,1], value=level)
contour_line.append(cl)
contour_lines.append(contour_line)
return contour_lines
def extract_contour_intervals(self, levels):
......@@ -586,14 +638,15 @@ class MeshGridData(GridData):
list with instances of :class:`MultiPolygonData`
"""
import numpy.ma as ma
import matplotlib._contour as _contour
from matplotlib._contour import QuadContourGenerator
values = ma.asarray(self.values)
contour_engine = _contour.QuadContourGenerator(self.lons, self.lats,
values.filled(), None, False, 0)
contour_engine = QuadContourGenerator(self.center_lons, self.center_lats,
values.filled(), None, False, 0)
contour_polygons = []
for lower_level, upper_level in zip(levels[:-1], levels[1:]):
segs, path_codes = contour_engine.create_filled_contour(lower_level, upper_level)
segs, path_codes = contour_engine.create_filled_contour(lower_level,
upper_level)
contour_mpg = MultiPolygonData([], [])
for i in range(len(segs)):
seg = segs[i]
......@@ -771,6 +824,14 @@ class GdalRasterData(MeshGridData):
function to apply to stored grid values.
Only applies if :param:`band_nr` is not null
(default: None)
Note:
Some differences with respect to :class:`MeshGridData`:
- Can be projected instead of lon/lat
- extent and resolution may be different from native extent and resolution
by setting a bounding box and down_sampling
- Native grid bounds and resolution in _x0, _x1, _y0, _y1, _dx, _dy
- Actual grid bounds and resolution in x0, x1, y0, y1, dx, dy
"""
def __init__(self, filespec, band_nr=1, down_sampling=1., nodata_value=None,
unit="", bbox=[], region=[], value_conversion=None):
......@@ -952,6 +1013,20 @@ class GdalRasterData(MeshGridData):
return int(round(abs((self.y0 - self.y1) / self.dy))) + 1
def get_bbox_from_region(self, region, margin_fraction=1./20):
"""
Construct bounding box from geographic region spec
:param region:
(min_lon, max_lon, min_lat, max_lat) tuple
:param margin_fraction:
float, width of margin as fraction of lon/lat range to
add around region to account for misalignment with
native coordinates
(default: 1./20)
:return:
list of floats: [llx, lly, urx, ury] in native coordinates
"""
from mapping.geotools.coordtrans import transform_coordinates
srs = self.srs
......@@ -1255,22 +1330,67 @@ class GdalRasterData(MeshGridData):
return out_data
def interpolate_grid(self, xout, yout, srs=None, checkbounds=False, masked=True, order=1):
values = self.interpolate(xout, yout, srs=srs, checkbounds=checkbounds, masked=masked, order=order)
def interpolate_grid(self, xout, yout, srs=None, checkbounds=False,
masked=True, order=1):
"""
Interpolate grid to a new grid.
Similar to :meth:`interpolate`, except that interpolation locations
should define a grid
:param:`xout`:
2-D array, X coordinates in native SRS or :param:`srs`
:param yout:
2-D array, Y coordinates in native SRS or :param:`srs`
:param srs:
:param checkbounds:
:param masked:
:param order:
see :meth:`interpolate`
:return:
instance of :class:`MeshGridData`
"""
values = self.interpolate(xout, yout, srs=srs, checkbounds=checkbounds,
masked=masked, order=order)
return MeshGridData(xout, yout, values)
def cross_section(self, xy0, xy1, num_points=100):
def cross_section(self, xy0, xy1, num_points=100, order=1):
"""
Native coordinates!
Extract cross section through grid.
:param xy0:
(X, Y) tuple defining start point in native grid coordinates
:param xy1:
(X, Y) tuple defining end point in native grid coordinates
:param num_points:
int, number of points in cross section
:param order:
see :meth:`interpolate`
:return:
1D array, interpolated grid values along cross section
"""
x0, y0 = xy0
x1, y1 = xy1
X = np.linspace(x0, x1, num_points)
Y = np.linspace(y0, y1, num_points)
## TODO: distances in km (take into account lon-lat rasters)
return self.interpolate(X, Y)
return self.interpolate(X, Y, order=order)
def warp_to_map(self, map, checkbounds=False, masked=True, order=1):
"""
Warp grid to map
:param map:
instance of :class:`LayeredBasemap`
:param checkbounds:
:param masked:
:param order:
see :meth:`interpolate`
:return:
2D array of warped grid coordinates
"""
from mapping.geotools.coordtrans import transform_mesh_coordinates
nx = int(round(map.figsize[0] * float(map.dpi)))
......@@ -1293,21 +1413,26 @@ class GdalRasterData(MeshGridData):
"""
from osgeo import gdal
## Open output format driver, see gdal_translate --formats for list
driver = gdal.GetDriverByName(format)
if driver and driver.GetMetadata().get(gdal.DCAP_CREATE) == 'YES':
## Open existing dataset
src_ds = gdal.Open(self.filespec, gdal.GA_ReadOnly)
if (self.down_sampling != 1
or self.get_bbox() != self.get_native_bbox()):
return self.to_mesh_grid().export_gdal(format, out_filespec)
else:
## Open output format driver, see gdal_translate --formats for list
driver = gdal.GetDriverByName(format)
if driver and driver.GetMetadata().get(gdal.DCAP_CREATE) == 'YES':
## Open existing dataset
src_ds = gdal.Open(self.filespec, gdal.GA_ReadOnly)
## Output to new format
dst_ds = driver.CreateCopy(out_filespec, src_ds, 0)
## Output to new format
dst_ds = driver.CreateCopy(out_filespec, src_ds, 0)
## Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
## Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
else:
print("Driver %s does not support CreateCopy() method" % format)
else:
print("Driver %s does not support CreateCopy() method" % format)
class MeshGridVectorData(BasemapData):
......
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