Interpolate dataset from TIN to regular mesh



In [ ]:

    
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter(action = "ignore", category = FutureWarning)

import os
import pandas as pd
import numpy as np

# Import badlands grid generation toolbox
import morphoGrid as morph
import pybadlands_companion.hydroGrid as hydr
from scipy.interpolate import RegularGridInterpolator

# display plots in SVG format
%config InlineBackend.figure_format = 'svg' 

%matplotlib inline

Interpolation from TIN

First we define the function buildRegularGrid that uses pyBadlands companion to remesh the triangular irregular network (TIN) over a regular grid.



In [ ]:

    
def buildRegularGrid(folder=None,step=None,outputname=None, basinXY=None):

    hydro = hydr.hydroGrid(folder=folder+'/h5', ncpus=1, \
                          ptXY = basinXY)
    hydro.getCatchment(timestep=step)
    
    morpho = morph.morphoGrid(folder=folder+'/h5', \
                          ncpus=1, dx=100)
    basinID = hydro.basinID
    morpho.loadHDF5(timestep=step,basinIDs=hydro.Basin, XYf=hydro.XY)
    
    df = pd.DataFrame({'X':morpho.x.flatten(),'Y':morpho.y.flatten(),
                        'Z':morpho.z.flatten(),'Basin':morpho.catchment.flatten()})   
    
    df.to_csv(outputname, columns=['X', 'Y', 'Z', 'Basin'], sep=',', index=False , 
              header=1)
    
    return basinID

Resolution of regular grid

One may want to change the resolution of the grid using the function regridDEM either to coarsen and refine the created grid.



In [ ]:

    
def regridDEM(res, disp, inDEM, outDEM):
    """
    Convert the initial regular DEM to the requested resolution.
    
    Parameters
    ----------
    variable: res
        Resampling resolution in metres.
    
    variable: disp
        Edge cell number to remove from the grid.
        
    variable: inDEM
        Name of the CSV topographic file to regrid.
        
    variable: outDEM
        Name of the new CSV topographic file.
    """

    xyz = pd.read_csv(str(inDEM), sep=',', engine='c', header=0, na_filter=False, \
                               dtype=np.float, low_memory=False)
    X = xyz.values[:,0]
    dx = X[1] - X[0]
    if res == dx:
        print 'Data spacing: ',dx
        print 'Requested spacing: ',res
        raise ValueError('The requested resolution is equal to the existing one.')
    Y = xyz.values[:,1]
    
    nx = int((X.max()-X.min())/dx + 1)
    ny = int((Y.max()-Y.min())/dx + 1)
    if nx*ny != len(X):
        raise ValueError('Check your input file the size of the grid is not right.')

    Z = np.reshape(xyz.values[:,2],(nx, ny),order='F')
    xgrid = np.arange(X.min(),X.max()+dx,dx)
    ygrid = np.arange(Y.min(),Y.max()+dx,dx)

    
    B = np.reshape(xyz.values[:,3],(nx, ny),order='F')
    RGI_function = RegularGridInterpolator((xgrid, ygrid), Z)
    RGI_function2 = RegularGridInterpolator((xgrid, ygrid), B)

    ngridX = np.arange(X.min(),X.max(),res)
    ngridY = np.arange(Y.min(),Y.max(),res)
    xi, yi = np.meshgrid(ngridX, ngridY)
    nnx = xi.shape[0]
    nny = xi.shape[1]

    zi = RGI_function((xi[disp:nnx-disp,disp:nny-disp].flatten(),
                       yi[disp:nnx-disp,disp:nny-disp].flatten()))
    
    bi = RGI_function2((xi[disp:nnx-disp,disp:nny-disp].flatten(),
                       yi[disp:nnx-disp,disp:nny-disp].flatten())).astype(int)
    
    df = pd.DataFrame({'X':xi[disp:nnx-disp,disp:nny-disp].flatten(),
                       'Y':yi[disp:nnx-disp,disp:nny-disp].flatten(),
                       'Z':zi.flatten(),'Basin':bi.flatten()})
    df.to_csv(str(outDEM),columns=['X', 'Y', 'Z','Basin'], sep=',', index=False)
    
    return

Exporting grids temporally

We first define the time steps number for which we wish to perform the regridding:



In [ ]:

    
tsteps = np.linspace(5,100,20).astype(int)

We then need to specify:

the regular grid resolution (dx),
the edge cell number to be removed (disps),
the pyBadlands output folder name containing the TIN files (foldername),
as well as an optional point coordinates (bXY) for which the corresponding basin index over time is going to be recorded.



In [ ]:

    
dx = 200.
disps = 5
bXY = [54331.5,4706.67]
foldername = '/workspace/volume/mountain/outputKeDecay'
outputfolder = '/workspace/volume/biodiversity/outputKeDecay'

We then loop over the requested time steps and perform the remeshing:



In [ ]:

    
basinID = []
for k in range(len(tsteps)):
    steps = tsteps[k]
    print '+ Performing regridding for step: ',steps
    output = 'tmpdata.csv'
    outputBio = foldername+str(steps)+'grid.csv'
    bID = buildRegularGrid(folder=foldername,step=steps,outputname=output,basinXY=bXY)
    basinID.append(bID)
    regridDEM(res = dx, disp = disps,inDEM = output, outDEM=outputBio)
    os.remove(output)
catchmentID = np.asarray(basinID,int)[:,0]

Catchment ID over time corresponding to the requested point:



In [ ]:

    
print catchmentID



In [ ]: