from __future__ import absolute_import, division, print_function, \
unicode_literals
import numpy as np
import pyflann
from scipy import spatial
import timeit
from rasterio.features import rasterize
from affine import Affine
import shapely.geometry
import shapely.ops
from functools import partial
from geometric_features.plot import subdivide_geom
from mpas_tools.mesh.creation.util import lonlat2xyz
[docs]def signed_distance_from_geojson(fc, lon_grd, lat_grd, earth_radius,
max_length=None, epsilon=1e-10):
"""
Get the distance for each point on a lon/lat grid from the closest point
on the boundary of the geojson regions.
Parameters
----------
fc : geometrics_features.FeatureCollection
The regions to be rasterized
lon_grd : numpy.ndarray
A 1D array of longitude values
lat_grd : numpy.ndarray
A 1D array of latitude values
earth_radius : float
Earth radius in meters
max_length : float, optional
The maximum distance (in degrees) between points on the boundary of the
geojson region. If the boundary is too coarse, it will be subdivided.
epsilon : float, optional
A small amount to expand each shape by to make sure lon/lat points on
the domain boundary are within shapes that touch the domain boundary.
Returns
-------
signed_distance : numpy.ndarray
A 2D field of distances (negative inside the region, positive outside)
to the shape boundary
"""
shapes = _subdivide_shapes(fc, max_length)
distance = distance_from_geojson(fc, lon_grd, lat_grd, earth_radius,
nn_search='flann', max_length=max_length,
shapes=shapes)
mask = mask_from_geojson(fc, lon_grd, lat_grd, shapes=shapes,
epsilon=epsilon)
signed_distance = (-2.0 * mask + 1.0) * distance
return signed_distance
[docs]def mask_from_geojson(fc, lon_grd, lat_grd, max_length=None, shapes=None,
epsilon=1e-10):
"""
Make a rasterized mask on a lon/lat grid from shapes (geojson multipolygon
data).
Parameters
----------
fc : geometrics_features.FeatureCollection
The regions to be rasterized
lon_grd : numpy.ndarray
A 1D array of longitude values
lat_grd : numpy.ndarray
A 1D array of latitude values
max_length : float, optional
The maximum distance (in degrees) between points on the boundary of the
geojson region. If the boundary is too coarse, it will be subdivided.
shapes : list of shapely.geometry, optional
A list of shapes that have already been extracted from fc and possibly
subdivided
epsilon : float, optional
A small amount to expand each shape by to make sure lon/lat points on
the domain boundary are within shapes that touch the domain boundary.
Returns
-------
mask : numpy.ndarray
A 2D mask with the shapes rasterized (0.0 outside, 1.0 inside)
"""
print("Mask from geojson")
print("-----------------")
if shapes is None:
shapes = _subdivide_shapes(fc, max_length)
nlon = len(lon_grd)
nlat = len(lat_grd)
uniform = (_is_uniform(lon_grd, epsilon=epsilon) and
_is_uniform(lat_grd, epsilon=epsilon))
if uniform:
dlon = (lon_grd[-1] - lon_grd[0])/(nlon-1)
dlat = (lat_grd[-1] - lat_grd[0])/(nlat-1)
# raserio works with pixels, and we want the lon/lat points to be at the
# centers of the pixels so the origin (the lower left of the first pixel) is
# half a pixel offset from the first lon/lat point
transform = Affine(dlon, 0., lon_grd[0] - 0.5*dlon,
0., dlat, lat_grd[0] - 0.5*dlat)
else:
shapes = _shapes_to_pixel_cooords(lon_grd, lat_grd, shapes)
transform = Affine(1., 0., 0.,
0., 1., 0.)
raster_shapes = []
for shape in shapes:
# expand a bit to make sure we hit the edges of the domain
shape = shape.buffer(epsilon)
raster_shapes.append((shapely.geometry.mapping(shape), 1.0))
mask = rasterize(raster_shapes, out_shape=(nlat, nlon), transform=transform)
return mask
[docs]def distance_from_geojson(fc, lon_grd, lat_grd, earth_radius, nn_search='flann',
max_length=None, shapes=None):
# {{{
"""
Get the distance for each point on a lon/lat grid from the closest point
on the boundary of the geojson regions.
Parameters
----------
fc : geometrics_features.FeatureCollection
The regions to be rasterized
lon_grd : numpy.ndarray
A 1D array of longitude values
lat_grd : numpy.ndarray
A 1D array of latitude values
earth_radius : float
Earth radius in meters
nn_search: {'kdtree', 'flann'}, optional
The method used to find the nearest point on the shape boundary
max_length : float, optional
The maximum distance (in degrees) between points on the boundary of the
geojson region. If the boundary is too coarse, it will be subdivided.
shapes : list of shapely.geometry, optional
A list of shapes that have already been extracted from fc and possibly
subdivided
Returns
-------
distance : numpy.ndarray
A 2D field of distances to the shape boundary
"""
print("Distance from geojson")
print("---------------------")
if shapes is None:
shapes = _subdivide_shapes(fc, max_length)
print(" Finding region boundaries")
boundary_lon = []
boundary_lat = []
for shape in shapes:
# get the boundary of each shape
x, y = shape.boundary.coords.xy
boundary_lon.extend(x)
boundary_lat.extend(y)
boundary_lon = np.array(boundary_lon)
boundary_lat = np.array(boundary_lat)
# Remove point at +/- 180 lon and +/- 90 lat because these are "fake".
# Need a little buffer (0.01 degrees) to avoid misses due to rounding.
mask = np.logical_not(np.logical_or(
np.logical_or(boundary_lon <= -179.99, boundary_lon >= 179.99),
np.logical_or(boundary_lat <= -89.99, boundary_lat >= 89.99)))
boundary_lon = boundary_lon[mask]
boundary_lat = boundary_lat[mask]
print(" Mean boundary latitude: {0:.2f}".format(np.mean(boundary_lat)))
# Convert coastline points to x,y,z and create kd-tree
npoints = len(boundary_lon)
boundary_xyz = np.zeros((npoints, 3))
boundary_xyz[:, 0], boundary_xyz[:, 1], boundary_xyz[:, 2] = \
lonlat2xyz(boundary_lon, boundary_lat, earth_radius)
flann = None
tree = None
if nn_search == "kdtree":
tree = spatial.KDTree(boundary_xyz)
elif nn_search == "flann":
flann = pyflann.FLANN()
flann.build_index(
boundary_xyz,
algorithm='kdtree',
target_precision=1.0,
random_seed=0)
else:
raise ValueError('Bad nn_search: expected kdtree or flann, got '
'{}'.format(nn_search))
# Convert backgound grid coordinates to x,y,z and put in a nx_grd x 3 array
# for kd-tree query
Lon_grd, Lat_grd = np.meshgrid(lon_grd, lat_grd)
X_grd, Y_grd, Z_grd = lonlat2xyz(Lon_grd, Lat_grd, earth_radius)
pts = np.vstack([X_grd.ravel(), Y_grd.ravel(), Z_grd.ravel()]).T
# Find distances of background grid coordinates to the coast
print(" Finding distance")
start = timeit.default_timer()
distance = None
if nn_search == "kdtree":
distance, _ = tree.query(pts)
elif nn_search == "flann":
_, distance = flann.nn_index(pts, checks=2000, random_seed=0)
distance = np.sqrt(distance)
end = timeit.default_timer()
print(" Done")
print(" {0:.0f} seconds".format(end-start))
# Make distance array that corresponds with cell_width array
distance = np.reshape(distance, Lon_grd.shape)
return distance
def _subdivide_shapes(fc, max_length):
shapes = []
for feature in fc.features:
# get the boundary of each shape
shape = shapely.geometry.shape(feature['geometry'])
if max_length is not None:
# subdivide the shape if it's too coarse
geom_type = shape.geom_type
shape = subdivide_geom(shape, geom_type, max_length)
shapes.append(shape)
return shapes
def _is_uniform(vector, epsilon=1e-10):
d = vector[1:] - vector[0:-1]
diff = d - np.mean(d)
return np.all(np.abs(diff) < epsilon)
def _shapes_to_pixel_cooords(lon, lat, shapes):
intx = partial(_interpx, lon)
inty = partial(_interpy, lat)
new_shapes = []
for shape in shapes:
shape = shapely.ops.transform(intx, shape)
shape = shapely.ops.transform(inty, shape)
new_shapes.append(shape)
return new_shapes
def _interpx(lon, x, y):
nlon = len(lon)
lon_pixels = np.arange(nlon, dtype=float)
x = np.interp(x, lon, lon_pixels)
return x, y
def _interpy(lat, x, y):
nlat = len(lat)
lat_pixels = np.arange(nlat, dtype=float)
y = np.interp(y, lat, lat_pixels)
return x, y