add sands case

mepelssystem/preprocessing/croco/preproc_croco_sands.m

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Build a CROCO RIP configuration.
+%  Create a grid file for the Rip Current experiment
+%  (for the parent and the child grid).
+%  This file is based on make_vortex.m
+%
+%  Further Information:
+%  http://www.croco-ocean.org
+%
+%  This file is part of CROCOTOOLS
+%
+%  CROCOTOOLS is free software; you can redistribute it and/or modify
+%  it under the terms of the GNU General Public License as published
+%  by the Free Software Foundation; either version 2 of the License,
+%  or (at your option) any later version.
+%
+%  CROCOTOOLS is distributed in the hope that it will be useful, but
+%  WITHOUT ANY WARRANTY; without even the implied warranty of
+%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%  GNU General Public License for more details.
+%
+%  You should have received a copy of the GNU General Public License
+%  along with this program; if not, write to the Free Software
+%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+%  MA  02111-1307  USA
+%
+%  Sandbar configuration
+%
+%  Julien CHAUCHAT, LEGI 2019
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+clear all
+close all
+
+
+%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Title
+%
+ROMS_title='SANDS';
+
+%
+% choose configuration between 'SANDS1', 'SANDS2'
+%
+mycase = "SANDS2"
+
+% choose input data between 'hydro (i.e. intermediate bed)', 'morpho (i.e. initial bed)'
+simulation = "morpho";
+%
+% choose input data between 'Polynomials', 'Measured'
+%
+inputDATA = "Measured"
+filterDATA = 1
+%
+%  Names
+%
+baseName = 'sandbar';
+%
+% grid size
+% (also modify Length_ETA accordingly)
+resolution = 'dx1m';
+gridSize = 1.;
+
+%
+% number of vertical layers
+%
+Nlayer=20;
+%
+% Variable roughness
+%
+zobOnshore = 2e-3;
+zobOffshore = 2e-3;
+xtransition = 150;
+delta = 2;
+
+% zobOnshore = 7e-3;
+% zobOffshore = 7e-4;
+% xtransition = 150;
+% delta = 2;
+%
+% max depth
+%
+depth=2.47;
+%
+% number of sediment class (not implemented yet)
+%
+Ns=1;
+%
+% number of bed layers
+%
+Sb=1;
+
+writeINI = 0;
+
+if mycase == 'SANDS1'
+    basePath= strcat('./../../beaches/sands/',simulation,'/sands1/');
+    
+    %/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip
+    %load('/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip/explip1b.mat');
+    %http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/explip1b.mat');
+    %s = lip1b;
+    
+    if simulation == 'hydro'
+        %nc=netcdf('http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/bathy_lip1b_inter.nc','nowrite');
+        nc = netcdf('/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip/bathy_lip1b_inter.nc','nowrite');
+
+        Xraw = nc{'X'}(:,:);
+        Yraw = nc{'Y'}(:,:);
+        Zraw = nc{'h'}(:,:);
+        
+        close(nc)
+        bathy = -Zraw;
+        
+    elseif simulation == 'morpho'
+        %nc=netcdf('http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/bathy_lip1b_init.nc','nowrite');
+        nc = netcdf('../../inputdata/bathy/sands/sands1_nonFilter_bathy.nc','nowrite');
+        
+       Xraw = nc{'X'}(:,:);
+        Yraw = nc{'Y'}(:,:);
+        Zraw = nc{'h'}(:,:);
+        
+        close(nc)
+        bathy = -Zraw;
+    end
+    
+else
+    basePath= strcat('./../../beaches/sands/',simulation,'/sands2/');
+    
+    %/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip
+    %load('/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip/explip1c.mat');
+    %http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/explip1c.mat');
+    %s = lip1c;
+    
+    if simulation == 'hydro'
+        %nc=netcdf('http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/bathy_lip1c_inter.nc','nowrite');
+        nc = netcdf('/fsnet/project/meige/2019/19MEPELS/OPENDATA_MEPELS/lip/bathy_lip1c_inter.nc','nowrite');
+
+        Xraw = nc{'X'}(:,:);
+        Yraw = nc{'Y'}(:,:);
+        Zraw = nc{'h'}(:,:);
+        
+        close(nc)
+        bathy = -Zraw;
+        
+    elseif simulation == 'morpho'
+        %nc=netcdf('http://servdap.legi.grenoble-inp.fr/opendap/meige/19MEPELS/lip/bathy_lip1c_init.nc','nowrite');
+        nc = netcdf('../../inputdata/bathy/sands/sands2_nonFilter_bathy.nc','nowrite');
+
+        Xraw = nc{'X'}(:,:);
+        Yraw = nc{'Y'}(:,:);
+        Zraw = nc{'h'}(:,:);
+        
+        close(nc)
+        bathy = -Zraw;
+    end
+end
+xbathy = Xraw;
+
+if (filterDATA==1)
+    dx=xbathy(1,2)-xbathy(1,1);
+    d1 = designfilt('lowpassiir','FilterOrder',12, ...
+        'HalfPowerFrequency',1./(2./dx),'DesignMethod','butter');
+    for j=1:size(bathy,1)
+        bathyf(j,:) = filtfilt(d1,bathy(j,:));
+    end
+else
+    bathyf = bathy;
+end
+%
+% Filenames
+%
+grdname = char(strcat(basePath,baseName,mycase,resolution,'_GRD.nc'));
+ininame = char(strcat(basePath,baseName,mycase,resolution,'_INI.nc'));
+
+%
+% min depth
+%
+hc=-20;
+
+zeta0=0;
+zetaMax=1;
+
+%
+% Parameters for the Grid
+%
+Length_XI =57;
+Length_ETA=gridSize;
+Xmin= -48.5;
+Ymin= 0;
+
+Nx= floor(Length_XI./gridSize)+2;
+Ny= floor(Length_ETA./gridSize)+2;
+
+dx=gridSize;
+dy=gridSize;
+
+Xmax= Xmin-dx/2 + dx*Nx + dx/2;
+Ymax= Ymin-dy/2 + dy*Ny + dy/2;
+
+Zmin= 0;
+Zmax= depth;
+%
+%
+% Horzontal Grid
+%
+%!----------------------------------------------------------------------
+%! Setup rectangulag grid:  coordinates (XI,ETA) at PSI- and RHO-points.
+%!----------------------------------------------------------------------
+for j=1:Ny
+    for i=1:Nx
+        xr(j,i) = Xmin + dx*((i-1-1)+0.5);
+        yr(j,i) = Ymin + dy*((j-1-1)+0.5);
+    end
+end
+%%min(min(xr));
+%%max(max(xr));
+
+X=xr;
+Y=yr;
+h=zeros(Ny,Nx);
+
+
+for j=1:Ny
+    for i=1:Nx
+        xx=xr(j,i);
+        if (mycase=="SANDS1")
+            if (inputDATA=="Measured")
+                dum = interp1(xbathy(1,:),bathyf(1,:),xx,'spline','extrap');
+                if (xx<min(xbathy(1,:)))
+                    dum=-depth;
+                end
+            end
+        elseif (mycase=="SANDS2")
+            if (inputDATA=="Measured")
+                dum = interp1(xbathy(1,:),bathyf(1,:),xx,'spline','extrap');
+                if (xx<min(xbathy(1,:)))
+                    dum=-depth;
+                end
+            end
+        end
+        h(j,i) = min([-dum, depth]);
+    end
+end
+%
+% Plot
+%
+
+fig1=figure;
+plot(xr(1,:),-h(1,:),'-g','LineWidth',1);
+if (inputDATA=="Measured")
+    hold on;
+    plot(xbathy(:,1),bathyf(:,1),'-k','LineWidth',1);
+    plot(xbathy(:,1),bathy(:,1),'-.r','LineWidth',1);
+end
+%axis([Xmin Xmax Zmin 770 -7 1 ]);
+xlabel('cross-shore (m)');
+ylabel('depth (m)');
+%axis([0,205,1.1, 4.2])
+legend('interpolated','filtered','raw','location','northwest')
+%%
+%%  create grid
+%%
+%     [Lonu,Lonv,Lonp]=rho2uvp(Lonr);
+%     [Latu,Latv,Latp]=rho2uvp(Latr);
+[xu,xv,xp]=rho2uvp(xr);
+[yu,yv,yp]=rho2uvp(yr);
+%
+% Create the grid file
+%
+disp(' ')
+disp(' Create the grid file...')
+[M,L]=size(yp);
+disp([' LLm = ',num2str(L-1)])
+disp([' MMm = ',num2str(M-1)])
+create_grid(L,M,grdname,ROMS_title)
+%
+% Fill the grid file
+%
+%     disp(' ')
+%     disp(' Fill the grid file...')
+%     nc=netcdf(grdname,'write');
+%     nc{'lat_u'}(:)=Latu;
+%     nc{'lon_u'}(:)=Lonu;
+%     nc{'lat_v'}(:)=Latv;
+%     nc{'lon_v'}(:)=Lonv;
+%     nc{'lat_rho'}(:)=Latr;
+%     nc{'lon_rho'}(:)=Lonr;
+%     nc{'lat_psi'}(:)=Latp;
+%     nc{'lon_psi'}(:)=Lonp;
+%     close(nc)
+%
+%  Compute the metrics
+%
+disp(' ')
+disp(' Compute the metrics...')
+[pm,pn,dndx,dmde]=get_metrics(grdname);
+pm=1./dx;
+pn=1./dy;
+dndx=0.;
+dmde=0.;
+
+xr=X;
+yr=Y;
+
+dx=1./pm;
+dy=1./pn;
+dxmax=max(max(dx/1000));
+dxmin=min(min(dx/1000));
+dymax=max(max(dy/1000));
+dymin=min(min(dy/1000));
+disp(' ')
+disp([' Min dx=',num2str(dxmin),' km - Max dx=',num2str(dxmax),' km'])
+disp([' Min dy=',num2str(dymin),' km - Max dy=',num2str(dymax),' km'])
+%
+%  Angle between XI-axis and the direction
+%  to the EAST at RHO-points [radians].
+%
+%angle = 0*angle ;
+rotat_angle = 0;
+%
+%  Coriolis parameter
+%
+%f=4*pi*sin(pi*Latr/180)/(24*3600);
+%f=ones(size(Latr))*1.04510e-4;
+f=ones(size(xr))*0.;
+%
+% Fill the grid file
+%
+disp(' ')
+disp(' Fill the grid file...')
+nc=netcdf(grdname,'write');
+nc{'pm'}(:)=pm;
+nc{'pn'}(:)=pn;
+nc{'dndx'}(:)=dndx;
+nc{'dmde'}(:)=dmde;
+nc{'x_u'}(:)=xu;
+nc{'y_u'}(:)=yu;
+nc{'x_v'}(:)=xv;
+nc{'y_v'}(:)=yv;
+nc{'x_rho'}(:)=xr;
+nc{'y_rho'}(:)=yr;
+nc{'x_psi'}(:)=xp;
+nc{'y_psi'}(:)=yp;
+nc{'angle'}(:)=rotat_angle;
+nc{'f'}(:)=f;
+nc{'h'}(:)=h;
+nc{'spherical'}(:)='F';
+%
+% Compute the mask
+%
+% mettre un facteur de secu
+maskr=h+zetaMax>0;
+maskr=process_mask(maskr);
+% test without mask
+%maskr(:,:)=1;
+[masku,maskv,maskp]=uvp_mask(maskr);
+%
+%  Write it down
+%
+nc{'h'}(:)=h;
+nc{'mask_u'}(:)=masku;
+nc{'mask_v'}(:)=maskv;
+nc{'mask_psi'}(:)=maskp;
+nc{'mask_rho'}(:)=maskr;
+
+%
+[Mp,Lp]=size(h);
+L=Lp-1;
+M=Mp-1;
+Np=Nlayer+1;
+
+%%%     Implementing non-uniform bed roughness %%%
+% zobt=7e-3;
+nc{'z0b'} = ncdouble('eta_rho','xi_rho') ;
+nc{'z0b'}.long_name = ncchar('bottom roughness');
+nc{'z0b'}.long_name = 'bottom roughnesss';
+nc{'z0b'}.units = ncchar('meter');
+nc{'z0b'}.units = 'meter';
+
+zob = zeros(Mp,Lp);
+% step function
+% zob2 = zeros(Mp,Lp);
+% zob2(:,:) = 1e-4;
+% zob2(find(xr(:,:)>140 & xr(:,:)<180)) = 7e-3;
+% hyperbolic tangent
+zob(:,:) = zobOffshore + (zobOnshore-zobOffshore)*0.5*(1+tanh((xr(:,:)-xtransition)/delta));
+
+nc{'z0b'}(:) = zob;
+
+figure
+plot(xr(1,:),zob(1,:),'--r')
+%hold on
+%plot(xr(1,:),zob2(1,:),'--r')
+xlabel('xr (m)')
+ylabel('zob (m)')
+
+%%%%%%
+%
+% nc{'z0b'} = ncdouble('eta_rho','xi_rho') ;
+% nc{'z0b'}.long_name = ncchar('bottom roughness');
+% nc{'z0b'}.long_name = 'bottom roughnesss';
+% nc{'z0b'}.units = ncchar('meter');
+% nc{'z0b'}.units = 'meter';
+% zob = zeros(Mp,Lp);
+% zob(:,:) = repmat(zobt,Mp,Lp);
+% nc{'z0b'}(:) = zob;
+%
+close(nc);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%                 create initial condition file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if (writeINI==1)
+    %
+    % Initialize the file
+    %
+    create_inifile(ininame,grdname,ROMS_title,...
+        0,0,hc,N,0,'clobber',1)
+    %
+    %  Read the grid file
+    %
+    nc=netcdf(grdname,'r');
+    h=nc{'h'}(:);
+    mask=nc{'mask_rho'}(:);
+    close(nc);
+    hmin=min(min(h(mask==1)));
+    
+    [Mp,Lp]=size(h);
+    L=Lp-1;
+    M=Mp-1;
+    Np=Nlayer+1;
+    %
+    %  open the initial file
+    %
+    nc = netcdf(ininame,'write');
+    %
+    % Write variables
+    %
+    time = 0.;%timeF(min(keepData));
+    nc{'spherical'}(:)='F';
+    nc{'tstart'}(:) =  time;
+    nc{'tend'}(:) =  time;
+    nc{'scrum_time'}(1) =  time*24*3600;
+    nc{'ocean_time'}(1) =  time*24*3600;
+    %
+    %     zobt=2e-2;
+    
+    nc{'zob'} = ncdouble('eta_rho','xi_rho') ;
+    nc{'zob'}.long_name = ncchar('bottom roughness');
+    nc{'zob'}.long_name = 'bottom roughnesss';
+    nc{'zob'}.units = ncchar('meter');
+    nc{'zob'}.units = 'meter';
+    zob = zeros(Mp,Lp);
+    zob(:,:) = repmat(zobt,Mp,Lp);
+    nc{'zob'}(:) = zob;
+    %
+    %     nc('s_b') = 1;
+    %
+    %     nc{'bed_thick'} = ncdouble('time','s_b','eta_rho','xi_rho') ;
+    %     %
+    %     nc{'bed_thick'}.long_name = ncchar('bed-layer-thickness');
+    %     nc{'bed_thick'}.long_name = 'bed layer thickness';
+    %     nc{'bed_thick'}.units = ncchar('meter');
+    %     nc{'bed_thick'}.units = 'meter';
+    %     bed_thick=zeros(1,1,Mp,Lp);
+    %     bed_thick(1,1,:,:) = repmat(2,Mp,Lp).*mask;
+    %     nc{'bed_thick'}(:) = bed_thick;
+    %
+    %
+    % Synchronize on disk
+    %
+    close(nc);
+end
+
+
+return

mepelssystem/preprocessing/mnt/configurations/sands_configuration.yaml

+# Truc Vert config file
+data_path:
+    local_Path: '/fsnet/project/meige/2019/19MEPELS/SANDS_WISE/bathy_sands1_init.nc'
+    
+# if you don't use a web link, leave the '' separate by a space else it could be a problem
+    web_path: ' '
+    extension: 'nc'
+    delimponct: "\t"
+
+    missingBathyPath: '/fsnet/project/meige/2019/19MEPELS/SANDS_WISE/bathy_sands1_init.nc'
+    delimPonctMissVal: "\t"
+
+    savePath_preproc: './../../inputdata/bathy/sands/'
+    
+
+parameter:
+    date_bathy: '1'
+# size of one side of one grid
+    gridsize: 0.5
+
+    colormaps: 'gist_earth'
+
+# to calculate the gradient bed elevation
+# 0 : no calcul / 1 : calcul
+    calGradH: 0
+
+# choice wath you want to save
+# 0 : just figures
+# 1 : just new bathy data
+# 2 : both
+    savePlotBathy: 2
+
+interpolation:
+# define the method used in the interpolation
+# for exemple : 'linear', 'quadratic', 'cubic', 'nearest'
+    interpMethod: 'linear'
+    
+# define the number of subsection on the bathymetry
+# the values define the number of subsection on one axis, so the real number is the square of your choice: ex : partition = 2 -> 4 subsection at all
+    partition: 1
+
+# you can save the interpolation data of raw bathymetry for an other calculate
+# 0 : no save after interpolation / 1 : save after interpolation 
+    savehraw: 1
+
+# in case you already done an interpolation with the same parameters before interpolation
+# 0 : not using backup / 1 : using backup
+    backupInterp: 0
+
+# path to the backup folder interpolation
+    folderBackupInterp: './backupInterp/'
+
+
+z0level:
+# mesure systeme of the z = 0 
+# zero hydrographique : put "hyd"
+# zero zero tide : put "tide"
+# zero IGN69 : put "ign"
+    zerolevel: "ign"
+# in case zerolevel is'nt zero hydrographique
+# give the difference (absolute) in meter to arrive at zero hydrographique
+    difflevel: 0
+
+
+missingValues:
+# method for replace mising values into bathymetry
+# 0 : no replacement / 1 : replacement with an other bathymetry / 2 : replacement with approximation
+    missVal: 2
+# if the study area miss more than a value, the bathymetry is automaticaly changed.
+# this paramter is involved just if you've choice if 1 for missVal
+    missmore: 100000
+# define the number of subsection on the bathymetry
+# the values define the number of subsection on one axis, so the real number is the square of your choice: ex : partition = 2 -> 4 subsection at all
+    partition: 2
+# defined the number of points used to made the approximation
+    ptsapprox: 4
+# define the method used in the interpolation
+# for exemple : 'linear', 'quadratic', 'cubic', 'nearest'
+    interpMethod: 'linear'
+
+changeStudyArea:
+# 0: enter the parameters directly in the code
+# 1: enter the parameters below
+    changeArea: 1
+    # point at the sea side
+    Xorig: -50.0
+    Yorig: 0.0
+    # cross-shore orientation
+    Lx: 58.0
+    # long-shore orientation
+    Ly: 0.5
+    origstudy: 'left'
+
+
+changeRotation:
+# defined the angle betwin the cross-shore and the north (the angle is the absolute value)
+    axisrotation: 0
+
+# if the bathymetry is not straight, to accord the study area "box" with the local landmark
+# 0 : inactivated / 1 : activated
+    rotArea: 0
+
+# give the rotation sense : to straighten the beach
+# 0 : trigo sense / 1 : anti trigo sense
+    trigosense: 0.0
+
+
+# if you want modify the coordinate system of the bathymetry
+projection:
+# the next 3 parameters are involved before interpolation
+# 0 : no transformation / 1 : transformation
+    before: 0
+# coordinate system bathymetry input
+    coorIn_B: "EPSG:2154"
+# final coordinate system
+    coorOut_B: "EPSG:2154"
+# EPSG: 
+## 2154 -> lambert93
+## 27571 -> lambert1
+## 27572 -> lambert2
+## 27573 -> lambert3
+## 27574 -> lambert4
+## 4326 -> WGS
+
+# the next 3 parameters are involved after interpolation and replacing values
+    after: 0
+    coorIn_A: "EPSG:2154"
+    coorOut_A: "EPSG:2154"
+
+# to change the geographic coordinate to local coordinate (beache)
+    changeCoordGeo: 0
+
+
+
+################################################################
+# input file for analyse and modif preprocessing mnt
+################################################################
+data_path_saved:
+    inputPath: "../../inputdata/bathy/sands/sands1_init.nc"
+
+
+#############################################################
+# analayse_preprocessing_mnt
+
+# conditions meanning on an a hard event onto 12h in the years 
+storm:
+    T12_Y: 16
+    H12_Y: 8.5
+# mean values on some months or also years 
+# not mandatory
+mean:
+    T_mean: 8
+    H_mean: 0.64
+# size of the sand grains
+sand_size:
+    D50: 250e-6
+
+sealevel:
+# to defined the sea level you've 3 choices: 
+# - IGN69: 0
+# - Middle-tide: 1
+# - zero-hydro: 2
+    defined_zero: 0
+# if choice different than zero hydraugraphique please give the 
+# diffence heigth between your choice and zero hydro 
+# if you don't know the delta between the h zero data and the zero-hydro
+# put : 'NaN'
+    diff_to_zero_hydro: 0
+# if you've the elevation for HTL and LTL fill the below variables
+# else put : 'NaN'   
+    HTL: 0
+    LTL: 0
+# approximation with the max height of the bathymetry 0 to 1
+    HTL_approx: 0.
+    tidalrange: 'nan'
+
+fit_dissipative:
+    sigma: 3
+    nbs_fit: 3
+
+#############################################################################
+# modif_bathy_preproc
+
+sandbar:
+    # possibility to create 0 / 1 or 2 bar
+    number_of_bar: 2
+    outerbar:
+        width: 150
+        height: 3
+        position: 810
+          # put 'NaN' if you want just a straight bar
+        lambda: NaN
+    innerbar:
+        width: 50
+        height: 0.5
+        position: 1190
+        lambda: NaN
+