# -*- coding: utf-8 -*-
"""
Mesh refinement functions for QMESH mesh gradation raster generation.
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
from shutil import copyfile
# Non-Standard Python Dependencies
import numpy as np
from scipy.interpolate import griddata
import shapefile
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
# Local Module Dependencies
from ocmw.core.graphics import gen_parula_cmap
from ocmw.dataman.dataReaders import ncread, ncwrite
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
[docs]
def writeRasterFile(rasterPath,baseRasterFile,rasterFile,raster):
"""
Write raster to netCDF file.
"""
copyfile(rasterPath+'/'+baseRasterFile,rasterPath+'/'+rasterFile)
ncwrite(rasterPath+'/'+rasterFile,'z',raster)
return
[docs]
def getRasterCoords(rasterPath,rasterFile):
"""
Get raster coordinates from netCDF file and generate gidded coordinates.
"""
xm = ncread(rasterPath+'/'+rasterFile,'x')
ym = ncread(rasterPath+'/'+rasterFile,'y')
[xq, yq] = np.meshgrid(xm,ym)
return [xm,ym,xq,yq]
[docs]
def normalise(v):
"""
Normalise a variable about its mean value.
"""
vnorm = v - np.nanmin(v)
vnorm = vnorm / np.nanmax(vnorm)
return vnorm
[docs]
def gradate(minval,maxval,var):
"""
Convert a raster into a gradation between minimum and maximum values.
"""
scl = var - np.nanmin(var[:])
delta = maxval-minval
grade = (delta/np.nanmax(scl[:]))*scl+minval;
return grade
[docs]
def gradation2meshDensity(gradation):
"""
Convert a gradation raster into an equivalent mesh density assuming equilateral triangles.
"""
mshDen = (4.0/np.sqrt(3))/np.square(gradation)
return mshDen
[docs]
def gradation2numElems(gradation,dx,dy):
"""
Convert a gradation raster into an equivalent number of element assuming equilateral triangles.
"""
dA = dx * dy
mshDen = gradation2meshDensity(gradation)
numElems = np.int64(np.sum(mshDen*dA))
return numElems
[docs]
def distance2coast(shpPath,shpFile,BndIDs,x,y):
"""
Calculate the distance for a set of points (x,y) from features in a coastline shapefile.
"""
sfile = shapefile.Reader(shpPath+'/'+shpFile)
nLines = sfile.numShapes
edgeX = None
edgeY = None
for ibnd in range(nLines):
feature = sfile.shapeRecords()[ibnd]
ID = sfile.records()[ibnd].as_dict()['PhysID']
if ID in BndIDs:
coords = np.asarray(feature.shape.points)
if edgeX is None:
edgeX = coords[:,0].copy()
edgeY = coords[:,1].copy()
else:
edgeX = np.concatenate((edgeX,coords[:,0]))
edgeY = np.concatenate((edgeY,coords[:,1]))
dist = np.empty(x.shape)
for pt in range(len(x)):
dx = edgeX - x[pt]
dy = edgeY - y[pt]
dr = np.sqrt(np.square(dx) + np.square(dy))
dist[pt] = np.min(dr)
return dist
[docs]
def rasterise(x,y,v,xq,yq):
"""
Convert a variable on an unstructured mesh into a raster.
"""
raster = griddata((x,y),v,(xq,yq),'linear')
return raster
[docs]
def vargradient(x,y,v,xq,yq):
"""
Calculate the spatial gradient of a mesh variable on a regular raster.
"""
xm = xq[0,:]
ym = yq[:,0]
dx = np.mean(np.diff(xm))
dy = np.mean(np.diff(ym))
vgrid = griddata((x,y),v,(xq,yq),'linear')
[dvdx,dvdy] = np.gradient(vgrid,dx,dy)
vgrad = np.sqrt(np.square(dvdx)+np.square(dvdy))
return vgrad
[docs]
def plotRaster(x,y,raster,minval,maxval):
"""
Plot as raster as a colourised image.
"""
cmap = gen_parula_cmap()
fig = plt.figure(figsize=(9,5))
ax = fig.gca()
divider = make_axes_locatable(ax)
ax_cb = divider.new_horizontal(size="4%", pad=0.1)
fig = ax.get_figure()
fig.add_axes(ax_cb)
im = ax.imshow(raster,origin='lower',cmap=cmap,extent=(np.min(x),np.max(x),np.min(y),np.max(y)),vmin=minval,vmax=maxval)
plt.colorbar(im, cax=ax_cb)
return fig,ax,ax_cb
#--------------------------------------------------------------------------
# MAIN PROCESS
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------