Commit 8d1a419e authored by Kris Vanneste's avatar Kris Vanneste
Browse files

Use absolute value of grid resolution to avoid problems when aligning bounding...

Use absolute value of grid resolution to avoid problems when aligning bounding box with the grid if resolution is negative in get_grid function.
Added 'constrain_bbox' argument to get_grid, extract_grid_data, extract_grid_profile and extract_grid_profiles functions.

git-svn-id: https://svn.seismo.oma.be/svn/seismo/mapping/layeredbasemap/trunk@7601 40b490c5-b4d9-47cb-8714-9bef99b524d5
parent 60d0685f
......@@ -21,7 +21,8 @@ __all__ = ['get_grid', 'get_grid_info', 'extract_grid_data',
'extract_grid_profile', 'extract_grid_profiles']
def get_grid(grid_spec, bbox=None, resolution=0, lonlat=False):
def get_grid(grid_spec, bbox=None, resolution=0, lonlat=False,
constrain_bbox=True):
"""
Fetch grid from file system or server, optionally applying a bounding
box or downsampling (for very large grids)
......@@ -51,7 +52,7 @@ def get_grid(grid_spec, bbox=None, resolution=0, lonlat=False):
grd = None
if not resolution:
wcs = WCSConnection(WCS_URL)
resolution = min(wcs.get_cell_size(layer_name))
resolution = np.abs(min(wcs.get_cell_size(layer_name)))
grd_bbox = wcs.get_bbox(layer_name)
else:
......@@ -68,7 +69,7 @@ def get_grid(grid_spec, bbox=None, resolution=0, lonlat=False):
grd_bbox = grd.get_bbox()
if not resolution:
resolution = min(grd.dx, grd.dy)
resolution = np.abs(min(grd.dx, grd.dy))
## Transform coordinates if necessary
if bbox:
......@@ -148,7 +149,7 @@ def get_grid_info(grid_spec):
def extract_grid_data(grid_spec, xdata, ydata, resolution=0, interp_order=1,
lonlat=False):
lonlat=False, constrain_bbox=True):
"""
Extract data from a potentially very large grid
......@@ -169,13 +170,20 @@ def extract_grid_data(grid_spec, xdata, ydata, resolution=0, interp_order=1,
:param lonlat:
bool, indicating whether or not X and Y positions are in degrees
(True) or in native grid coordinates (False)
:param constrain_bbox:
bool, whether or not to apply a bounding box to the grid before
extracting the data in order to improve efficiency
(default: True)
:return:
numpy float array, grid data
"""
xdata = np.asarray(xdata)
ydata = np.asarray(ydata)
bbox = (xdata.min(), ydata.min(), xdata.max(), ydata.max())
if constrain_bbox:
bbox = (xdata.min(), ydata.min(), xdata.max(), ydata.max())
else:
bbox = None
grd = get_grid(grid_spec, bbox=bbox, resolution=resolution, lonlat=lonlat)
......@@ -190,7 +198,7 @@ def extract_grid_data(grid_spec, xdata, ydata, resolution=0, interp_order=1,
def extract_grid_profile(grid_spec, xy0, xy1, num_points=None, interp_order=1,
lonlat=False):
lonlat=False, constrain_bbox=True):
"""
Extract data from a grid along a profile (cross-section)
......@@ -211,6 +219,10 @@ def extract_grid_profile(grid_spec, xy0, xy1, num_points=None, interp_order=1,
:param lonlat:
bool, indicating whether or not X and Y positions are in degrees
(True) or in native grid coordinates (False)
:param constrain_bbox:
bool, whether or not to apply a bounding box to the grid before
extracting the data in order to improve efficiency
(default: True)
:return:
(X, Y, distances, data) tuple of 1D arrays
......@@ -221,10 +233,12 @@ def extract_grid_profile(grid_spec, xy0, xy1, num_points=None, interp_order=1,
"""
xmin, xmax = sorted([xy0[0], xy1[0]])
ymin, ymax = sorted([xy0[1], xy1[1]])
bbox = (xmin, ymin, xmax, ymax)
if constrain_bbox:
bbox = (xmin, ymin, xmax, ymax)
else:
bbox = None
grd = get_grid(grid_spec, bbox=bbox, resolution=0, lonlat=lonlat)
xmin, ymin, xmax, ymax = bbox
xdata = np.array([xmin, xmax])
ydata = np.array([ymin, ymax])
......@@ -254,7 +268,7 @@ def extract_grid_profile(grid_spec, xy0, xy1, num_points=None, interp_order=1,
def extract_grid_profiles(grid_spec, xy0, xy1, num_points=None, interp_order=1,
lonlat=False):
lonlat=False, constrain_bbox=True):
"""
Similar to func:`extract_grid_profile`, but for multiple profiles at once
......@@ -274,7 +288,10 @@ def extract_grid_profiles(grid_spec, xy0, xy1, num_points=None, interp_order=1,
Y = np.array([pt[1] for pt in xy0 + xy1])
xmin, xmax = X.min(), X.max()
ymin, ymax = Y.min(), Y.max()
bbox = (xmin, ymin, xmax, ymax)
if constrain_bbox:
bbox = (xmin, ymin, xmax, ymax)
else:
bbox = None
grd = get_grid(grid_spec, bbox=bbox, resolution=0, lonlat=lonlat)
......
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