Commit 10dff3f9 authored by Kris Vanneste's avatar Kris Vanneste
Browse files

Added grid.py containing ROB-specific grid functions.

Updated __init__.py.

git-svn-id: https://svn.seismo.oma.be/svn/seismo/mapping/layeredbasemap/trunk@7547 40b490c5-b4d9-47cb-8714-9bef99b524d5
parent 9411bcbd
......@@ -24,7 +24,14 @@ else:
## Import submodules
## geoserver
## grid
if not reloading:
from . import grid
else:
reload(grid)
from .grid import *
## geoserver (depends on grid)
if not reloading:
from . import geoserver
else:
......
# -*- coding: utf-8 -*-
"""
Custom functions related to grids on ROB-specific locations
(seismogis, geoserver, ...)
Created on Thu Feb 11 12:40:35 2021
@author: kris
"""
from __future__ import absolute_import, division, print_function, unicode_literals
import os
import numpy as np
from ..data_types import GdalRasterData, WCSData
__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):
"""
Fetch grid from file system or server, optionally applying a bounding
box or downsampling (for very large grids)
:param grid_spec:
str, full path to grid file or name of dataset known by seismogis
or name of layer on geoserver
:param bbox:
list or tuple of floats: (llx, lly, urx, ury), bounding box to clip grid
(default: None, will use entire grid)
:param resolution:
float, resolution of DEM in native units (for potential downsampling)
(default: 0, will use native grid resolution)
:param lonlat:
bool, indicating whether or not :param:`bbox`is in degrees
(True) or in native grid coordinates (False)
(default: False)
:return:
instance of :class:`GdalRasterData` or :class:`WCSData`
"""
if ':' in grid_spec[2:]:
from ..data_types.grid import WCSConnection
from .geoserver import WCS_URL
layer_name = grid_spec
grd = None
if not resolution:
wcs = WCSConnection(WCS_URL)
resolution = min(wcs.get_cell_size(layer_name))
grd_bbox = wcs.get_bbox(layer_name)
else:
if os.path.exists(grid_spec):
grd_file = grid_spec
else:
from mapping.seismogis import SeismoGisCatalog as giscat
giscat.load_collections(verbose=False)
[ds] = giscat.find_datasets(grid_spec)
grd_file = ds.get_gis_filespec()
grd = GdalRasterData(grd_file)
grd_bbox = grd.get_bbox()
if not resolution:
resolution = min(grd.dx, grd.dy)
## Transform coordinates if necessary
if bbox:
xmin, ymin, xmax, ymax = bbox
xdata = np.array([xmin, xmax])
ydata = np.array([ymin, ymax])
if lonlat:
from mapping.geotools.coordtrans import (transform_array_coordinates, WGS84)
if grd:
grd_srs = grd.srs
else:
grd_srs = wcs.get_srs(layer_name)
xdata, ydata = transform_array_coordinates(WGS84, grd_srs, xdata, ydata)
## Determine bounding box
xmin = np.floor(xdata.min() / resolution) * resolution
xmax = np.ceil(xdata.max() / resolution) * resolution
ymin = np.floor(ydata.min() / resolution) * resolution
ymax = np.ceil(ydata.max() / resolution) * resolution
xmin = max(xmin, grd_bbox[0])
xmax = min(xmax, grd_bbox[2])
ymin = max(ymin, grd_bbox[1])
ymax = min(ymax, grd_bbox[3])
bbox = (xmin, ymin, xmax, ymax)
if (xmin > xmax) or (ymin > ymax):
raise Exception('Invalid bounding box: ' ,bbox)
if grd:
if bbox:
grd.apply_bbox(bbox)
else:
grd = WCSData(WCS_URL, layer_name, resolution=resolution, bbox=bbox)
return grd
def get_grid_info(grid_spec):
"""
Get information about a grid without having to read the entire file
:param grid_spec:
str, full path to grid file or name of dataset known by seismogis
or name of layer on geoserver
:return:
dict, having 'srs', 'cell_size' and 'bbox' keys
"""
grd_info = {}
if ':' in grid_spec[2:]:
from ..data_types.grid import WCSConnection
from .geoserver import WCS_URL
layer_name = grid_spec
wcs = WCSConnection(WCS_URL)
grd_info['srs'] = wcs.get_srs(layer_name)
grd_info['cell_size'] = wcs.get_cell_size(layer_name)
grd_info['bbox'] = wcs.get_bbox(layer_name)
else:
if os.path.exists(grid_spec):
grd_file = grid_spec
else:
from mapping.seismogis import SeismoGisCatalog as giscat
giscat.load_collections(verbose=False)
[ds] = giscat.find_datasets(grid_spec)
grd_file = ds.get_gis_filespec()
grd = GdalRasterData(grd_file)
grd_info['srs'] = grd.srs
grd_info['cell_size'] = (abs(grd.dx), abs(grd.dy))
grd_info['bbox'] = grd.get_bbox()
return grd_info
def extract_grid_data(grid_spec, xdata, ydata, resolution=0, interp_order=1,
lonlat=False):
"""
Extract data from a potentially very large grid
:param grid_spec:
str, full path to grid file or name of dataset known by seismogis
or name of layer on geoserver
:param xdata:
list or array of floats, X positions in native coordinates of DEM
:param ydata:
list or array of floats, Y positions in native coordinates of DEM
:param resolution:
float, resolution of grid in native units
(default: 0, will use best resolution available)
:param interp_order:
int, type of interpolation, 0=nearest neighbor, 1=bilinear,
3=cubic spline
(default: 1)
:param lonlat:
bool, indicating whether or not X and Y positions are in degrees
(True) or in native grid coordinates (False)
:return:
numpy float array, grid data
"""
xdata = np.asarray(xdata)
ydata = np.asarray(ydata)
bbox = (xdata.min(), xdata.max(), ydata.min(), ydata.max)
grd = get_grid(grid_spec, bbox=bbox, resolution=resolution, lonlat=lonlat)
## Transform coordinates if necessary
if lonlat:
from mapping.geotools.coordtrans import (transform_array_coordinates, WGS84)
xdata, ydata = transform_array_coordinates(WGS84, grd.srs, xdata, ydata)
data = grd.interpolate(xdata, ydata, order=interp_order)
return data
def extract_grid_profile(grid_spec, xy0, xy1, num_points=None, interp_order=1,
lonlat=False):
"""
Extract data from a grid along a profile (cross-section)
:param grid_spec:
str, full path to grid file or name of dataset known by seismogis
or name of layer on geoserver
:param xy0:
(x, y) tuple, coordinates of starting point of profile
:param xy1:
(x, y) tuple, coordinates of ending point of profile
:param num_points:
int, number of points to discretize profile in
(default: None, will auto-determine based on native resolution of grid)
:param interp_order:
int, type of interpolation, 0=nearest neighbor, 1=bilinear,
3=cubic spline
(default: 1)
:param lonlat:
bool, indicating whether or not X and Y positions are in degrees
(True) or in native grid coordinates (False)
:return:
(X, Y, distances, data) tuple of 1D arrays
- X: X coordinates of profile points (native grid units or degrees)
- Y: Y coordinates of profile points (native grid units or degrees)
- distances: distances of profile points (native grid units)
- data: grid data along profile
"""
xmin, xmax = sorted([xy0[0], xy1[0]])
ymin, ymax = sorted([xy0[1], xy1[1]])
bbox = (xmin, ymin, xmax, ymax)
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])
## Transform coordinates if necessary
if lonlat:
from mapping.geotools.coordtrans import (transform_array_coordinates, WGS84)
xdata, ydata = transform_array_coordinates(WGS84, grd.srs, xdata, ydata)
xmin, xmax = xdata
ymin, ymax = ydata
if num_points is None:
resolution = min(abs(grd.dx), abs(grd.dy))
d = np.hypot((xmax-xmin), (ymax-ymin))
num_points = int(np.round(d / resolution) + 1)
xdata = np.linspace(xmin, xmax, num_points)
ydata = np.linspace(ymin, ymax, num_points)
distances = np.hypot((xdata-xmin), (ydata-ymin))
data = grd.interpolate(xdata, ydata, order=interp_order)
if lonlat:
xdata, ydata = transform_array_coordinates(grd.srs, WGS84, xdata, ydata)
return (xdata, ydata, distances, data)
def extract_grid_profiles(grid_spec, xy0, xy1, num_points=None, interp_order=1,
lonlat=False):
"""
Similar to func:`extract_grid_profile`, but for multiple profiles at once
:param xy0:
list of (x, y) tuples, profile starting points
:param xy1:
list of (x, y) tuples, profile ending points
:return:
list with (X, Y, distances, data) tuples of 1D arrays
"""
num_profiles = len(xy0)
assert len(xy1) == len(xy0)
## Determine bounding box
X = np.array([pt[0] for pt in xy0 + xy1])
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)
grd = get_grid(grid_spec, bbox=bbox, resolution=0, lonlat=lonlat)
if lonlat:
from mapping.geotools.coordtrans import (transform_array_coordinates, WGS84)
result = []
for i in range(num_profiles):
X = np.array([xmin, xmax])
Y = np.array([ymin, ymax])
## Transform coordinates if necessary
if lonlat:
X, Y = transform_array_coordinates(WGS84, grd.srs, X, Y)
xmin, xmax = X
ymin, ymax = Y
if num_points is None:
resolution = min(abs(grd.dx), abs(grd.dy))
d = np.hypot((xmax-xmin), (ymax-ymin))
num_points = int(np.round(d / resolution) + 1)
xdata = np.linspace(xmin, xmax, num_points)
ydata = np.linspace(ymin, ymax, num_points)
distances = np.hypot((xdata-xmin), (ydata-ymin))
data = grd.interpolate(xdata, ydata, order=interp_order)
if lonlat:
xdata, ydata = transform_array_coordinates(grd.srs, WGS84, xdata, ydata)
result.append((xdata, ydata, distances, data))
return result
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