import os
import numpy as np
from cem_gridtools.gridgen import gridgen
from cem_gridtools.grid_interp import init, on_point
from shapely.geometry import Point, Polygon
from shapely.geometry import MultiLineString, MultiPoint
import shapely
from vincenty import vincenty
# Create xy tuple from x and y arrays
def _xy(x,y):
return tuple([(x[i],y[i]) for i in range(x.shape[0])])
[docs]class CGrid(object):
""" CGrid constructor
Args:
x (array): x boundary values
y (array): y boundary values
beta (array): beta values (corner status)
*beta* points
* +1 for external corner or 90 degree clockwise turns
* -1 for internal corner or 90 degree anti-clockwise turns
* 0 all other boundary points
"""
def __init__(self,x,y,beta):
# Boundary values
self.x = x
self.y = y
# Corner mask
self.beta = beta
# grid output
self.gx = None
self.gy = None
# grid dimensions
self.ni = -1
self.nj = -1
# Double density grid
self.ddx = None
self.ddy = None
def __str__(self):
if self.ni > -1:
return "CGrid object : generated grid with ni={}, nj={}".format(self.ni, self.nj)
else:
return "CGrid object : no grid"
[docs] def gen_grid(self,ni,nj,ulidx=0):
""" Generates the **(ni+1) x (nj+1)** grid and calculates the full double-density points
Args:
ni (int): Number of required grid **cells** in the i-direction
nj (int): Number of required grid **cells** in the j-direction
ulidx (int, optional): The index of the *upper left* corner of the grid
"""
self.ni = ni+1 # one for grid points
self.nj = nj+1
# Call gridgen-c
self.gx, self.gy = gridgen(self.x, self.y, self.beta, self.ni, self.nj, ulidx)
# Helper function to generate double-denstity points
self.gen_ddgrid()
def gen_ddgrid(self):
gx = self.gx
gy = self.gy
# Define full shape
dshp = (2*self.gx.shape[0]-1, 2*self.gx.shape[1]-1)
ddx = np.zeros(dshp, dtype='d')
ddy = np.zeros(dshp, dtype='d')
# Fill the single density values
ddx[0::2, 0::2] = gx[:]
ddy[0::2, 0::2] = gy[:]
# Left face
ddx[0::2,1:-1:2] = 0.5*(gx[:,0:-1] + gx[:,1:])
ddy[0::2,1:-1:2] = 0.5*(gy[:,0:-1] + gy[:,1:])
# Back face
ddx[1:-1:2,0::2] = 0.5*(gx[0:-1,:] + gx[1:,:])
ddy[1:-1:2,0::2] = 0.5*(gy[0:-1,:] + gy[1:,:])
# Centres
ddx[1:-1:2, 1:-1:2] = 0.25*(gx[0:-1,0:-1] + gx[1:,1:] + gx[0:-1,1:] + gx[1:,0:-1])
ddy[1:-1:2, 1:-1:2] = 0.25*(gy[0:-1,0:-1] + gy[1:,1:] + gy[0:-1,1:] + gy[1:,0:-1])
# Set attributes
self.ddx = ddx
self.ddy = ddy
[docs] def plot(self, legend=False, dd=False, map=False):
""" Simple plotting routine for gridlines using Matplotlib
"""
import matplotlib.pyplot as plt
if map:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
proj = None
if map:
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=[8,6])
ax = fig.add_subplot(projection=proj)
if map:
coast = cfeature.GSHHSFeature()
ax.add_feature(coast, zorder=0)
ax.gridlines(draw_labels=True)
if self.gx is not None:
ax.plot(self.gx, self.gy, 'b-')
ax.plot(self.gx.T, self.gy.T, 'b-')
if dd:
ax.plot(self.ddx, self.ddy, 'go')
ax.plot(self.ddx.T, self.ddy.T, 'go')
# External points
I = self.beta>0
Le, = ax.plot(self.x[I], self.y[I], 'ks', label="External corner")
# Internal points
I = self.beta<0
Li, = ax.plot(self.x[I], self.y[I], 'k^', label="Internal corner")
# All others
I = self.beta==0
Lo, = ax.plot(self.x[I], self.y[I], 'k*', label='Other boundary')
if legend:
ax.legend(handles=[Le, Li, Lo], bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
ax.autoscale()
plt.show()
[docs] def get_grid_points(self):
""" Returns the 2D grid corner points (gx, gy)
"""
if self.gx is None:
raise Exception("Generate the grid first")
return self.gx, self.gy
[docs] def get_cell_centres(self):
""" Returns the 2D cell centres for (xc, yc) for the grid
"""
if self.ddx is None:
raise Exception("Generate the grid first")
xc = self.ddx[1:-1:2, 1:-1:2]
yc = self.ddy[1:-1:2, 1:-1:2]
return xc, yc
[docs] def export_shoc_grid(self, name):
"""Write SHOC grid to given file
Args:
name (str): output filename, overwrites file
"""
if self.ddx is None:
raise Exception("Generate the grid first")
N = self.ddx.size
with open(name, 'w') as f:
f.write('# X coordinates\n')
f.write('XCOORDS {}\n'.format(N))
np.savetxt(f, self.ddx.ravel(), fmt='%.7f')
f.write('\n')
f.write('# Y coordinates\n')
f.write('YCOORDS {}\n'.format(N))
np.savetxt(f, self.ddy.ravel(), fmt='%.7f')
[docs] def get_cell_centres_as_points(self):
""" Returns the cell centres as a **shapely** MultiPoint object
"""
if self.ddx is None:
raise Exception("Generate the grid first")
xc = self.ddx[1:-1:2, 1:-1:2]
yc = self.ddy[1:-1:2, 1:-1:2]
I = np.isfinite(xc.ravel())
return MultiPoint(_xy(xc.ravel()[I], yc.ravel()[I]))
[docs] def get_grid_as_lines(self):
""" Returns the grid as a **shapely** MultiLineString object
* This is still work in progress and the API is currently rapidly evolving ...
"""
if self.gx is None:
raise Exception("Generate the grid first")
gx, gy = self.gx, self.gy
lines = []
lines += [_xy(gx[i,:],gy[i,:]) for i in range(gx.shape[0])]
lines += [_xy(gx[:,j],gy[:,j]) for j in range(gx.shape[1])]
return MultiLineString(tuple(lines))
[docs] def get_grid_as_polygons(self):
""" Returns the grid as a list of **shapely** Polygon objects
A MultiPolygon would have invalid geometry as they cannot share an edge
"""
if self.gx is None:
raise Exception("Generate the grid first")
gx = self.gx
gy = self.gy
P = []
# C [row] - major order
for j in range(gx.shape[0]-1):
for i in range(gx.shape[1]-1):
# Each cell
X = gx[j:j+2,i:i+2]
Y = gy[j:j+2,i:i+2]
if np.any(np.isnan(X)):
continue
# Anti-clockwise points
idx = np.array([0,2,3,1,0])
X = X.ravel()[idx]
Y = Y.ravel()[idx]
poly = P.append(Polygon(_xy(X,Y)))
# Return as list
return P
[docs]class CGridRect(CGrid):
""" CGRidRect constructor
If *geo=True* then the resolution will be interpreted in "metres" from which the grid dimensions
are calcuated based on the geographic centroid of the grid. Otherwise, units are assummed to
be the same as the supplied input boundary
Args:
x (array): 1-d array of length 4
y (array): 1-d array of length 4
geo (bool, optional): whether this is geographic
"""
def __init__(self,x,y,geo=False):
super().__init__(x, y, np.ones(4, dtype='d'))
self.isgeo = geo
[docs] def gen_gridrect(self, resolution):
""" Generate uniform grid
Args:
resolution (float): Required grid resolution
"""
x,y = (self.x, self.y)
if len(x) != 4 and len(y) != 4:
raise Exception("x and y arrays must both be of length 4")
# Form Polygon to ensure anti clockwise orientation
pts = ((x[0],y[0]), (x[1],y[1]), (x[2],y[2]), (x[3],y[3]))
P = shapely.geometry.polygon.orient(Polygon(pts))
c = P.exterior.coords
# The first index is ul, the next down -j
if self.isgeo:
# Vincenty points are (lat,lon) so opposite of (x,y)!!
dist_i = vincenty((c[0][1], c[0][0]), (c[3][1], c[3][0])) * 1000 # convert to m
dist_j = vincenty((c[0][1], c[0][0]), (c[1][1], c[1][0])) * 1000
else:
dist_i = Point(c[0][0], c[0][1]).distance(Point(c[3][0], c[3][1]))
dist_j = Point(c[0][0], c[0][1]).distance(Point(c[1][0], c[1][1]))
self.ni = int(dist_i / resolution) + 1
self.nj = int(dist_j / resolution) + 1
# print('ni={}, nj={}'.format(self.ni, self.nj))
self.gx, self.gy = gridgen(self.x, self.y, self.beta, self.ni, self.nj, 0)
# Helper function to generate double-denstity points
self.gen_ddgrid()
# end CGridRect
[docs]class NN(object):
""" NN constructor
Args:
x (array): X points (1-d)
y (array): Y points (1-d)
z (array): Z (depth) points (1-d)
irule (string):
interpolation method to use
Choices are **nn_sibson** (default), **nn_non_sibson** and **linear**
"""
def __init__(self,x,y,z,irule="nn_sibson"):
self.h = init(x,y,z,irule)
self.V = None
[docs] def interp(self,xo,yo):
""" Perform interpolation on the given points
Args:
xo (float): x value(s) to interpolate onto
yo (float): y value(s) to interpolate onto
**xo,yo** may be scalars or 1-d arrays
Returns:
float: interpolated value(s)
"""
if isinstance(xo, np.ndarray):
# Check array
if len(xo.shape):
N = len(xo)
v = np.zeros(N, dtype='d')
for i in range(N):
v[i] = on_point(self.h, xo[i], yo[i])
self.V = v
return v
# Assume float scalar
v = on_point(self.h, xo, yo)
self.V = v
return v
[docs] def export_shoc_bathy(self, name, iV=None):
"""Write SHOC bathymetry data to given file
Args:
name (str): output filename, overwrites file
iV (array, optional): Input vector to export
"""
if self.V is None:
raise Exception("Perform interpolation first")
V = self.V
if iV is not(None):
if not(V.shape == iV.shape):
raise Exception("export_shoc_bathy: Incompatible shape of given vector")
V = iV
# Write to file
with open(name, 'w') as f:
f.write('\n')
f.write('# Bathymetry\n')
f.write('BATHY {}\n'.format(V.size))
np.savetxt(f, V.ravel(), fmt='%.2f')
f.write('\n')
[docs]class Coast(object):
""" Coast constructor
The type of data is inferred from the extension:
* .cst : ASCII X,Y columns with NaN breaks
* .shp : Shape file (*not implemented yet*)
Args:
fname (str): coast file name
"""
def __init__(self,fname):
self.pts = np.loadtxt(fname)
# Get all NaN indices
self.idx = np.flatnonzero(np.isnan(self.pts[:,0]))
[docs] def get_num_polygons(self):
""" Get the number of polygons found in the coasts file
Returns:
integer: number of polygons
"""
return self.idx.shape[0]
def get_polyons_as_generator(self):
from shapely.validation import explain_validity
""" Generator function for polygons in the coasts file
Example::
fname = '/home/farhan/work/data/coasts/aust_200m.cst'
C = Coast(fname)
polyFunc = C.get_polyons_as_generator()
poly = next(polyFunc) # Get first polygon
# To iterate over all polygons
while poly is not None:
# Input data with less that 2 points are initialised to empty
if not(poly.is_empty):
# Do something with 'poly'
poly = next(polyFunc) # continue onto next
"""
pts = self.pts
idx = self.idx
st = 0 # start index
for i in idx:
en = i # end index
if en-st > 3: # Make sure we have enough points
p = Polygon(pts[st:en,:])
if p.is_valid:
yield p
else:
print("Polygon using indices ({}:{}) is invalid:{}".format(st,en,explain_validity(p)))
if p.buffer(0).is_valid:
yield p.buffer(0)
yield Polygon()
else:
yield Polygon()
# print("st={} en={}".format(st,en))
st = en+1
# final
yield None
[docs] def plot(self, ax=None):
""" Plot all coastline polygons
Args:
ax (optional): Axis to plot on
"""
def get_xy_from_poly(poly):
x = []; y = []
for c in poly.exterior.coords:
x.append(c[0])
y.append(c[1])
return x,y
if ax is None:
fig, ax = plt.subplots()
pfunc = self.get_polyons_as_generator()
poly = next(pfunc)
while poly is not None:
if not(poly.is_empty):
if poly.geom_type == 'MultiPolygon':
for pp in poly.geoms:
x,y = get_xy_from_poly(pp)
ax.plot(x,y,'b')
else:
x,y = get_xy_from_poly(poly)
ax.plot(x,y,'b')
poly = next(pfunc)
plt.show()
# end cem_gridtools.py