Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 65bdba7a authored by Lucia Rinchiuso's avatar Lucia Rinchiuso
Browse files

Selection by ID, country and river

parents b43c409d 3378ab20
......@@ -9,16 +9,20 @@ FFLAG = -fPIC
FPY = f2py
#PDBG = --debug
#
EXT_SUFFIX := $(shell python3-config --extension-suffix)
#
all : routing_interface.so diagst.so
routing_interface : routing_interface$(EXT_SUFFIX)
diagst : diagst$(EXT_SUFFIX)
clean :
rm -f compilation.log *.o *.so *.mod routing_interface.so diagst.so *~
routing_interface.so : $(obj) routing_interface.f90
routing_interface$(EXT_SUFFIX) : $(obj) routing_interface.f90
rm -f compilation.log
$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=gfortran -m routing_interface $(obj) > compilation.log
diagst.so : diagst.f90
diagst$(EXT_SUFFIX) : diagst.f90
$(FPY) -c diagst.f90 $(PDBG) --fcompiler=gfortran -m diagst > compilation.log
defprec.o : defprec.f90
......
......@@ -165,7 +165,7 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
IF ( debug) WRITE(numout,*) "Memory Mgt findbasin : nbvmax, nb_htu, nbv = ", nbvmax, nb_htu, nbv
DO ib=1,nbpt
CALL routing_reg_findbasins(nb_htu, nbv, ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
CALL routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
& basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
& topoind_bx(ib,:,:), lshead_bx(ib,:,:), diaglalo, &
& nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), basin_sz(ib,:), basin_bxout(ib,:), &
......
......@@ -41,7 +41,7 @@ MODULE routing_reg
! values are then averaged.
! topo_option = 4 : As above but the average topoindex is computed for each river.
!
INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
INTEGER(i_std), SAVE :: topo_option = 3 !! Option to calculate topoindex
!
CONTAINS
!! ================================================================================================================================
......@@ -384,7 +384,7 @@ CONTAINS
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, &
SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, &
& nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, &
& basin_lshead, coast_pts, lontmp, lattmp)
!
......@@ -393,18 +393,19 @@ CONTAINS
!! INPUT VARIABLES
INTEGER(i_std), INTENT(in) :: nb_htu, nbv
INTEGER(i_std), INTENT(in) :: ib !!
INTEGER(i_std), INTENT(in) :: ijdimmax
INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless)
INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless)
REAL(r_std), INTENT(in) :: hierarchy(:,:) !!
REAL(r_std), INTENT(in) :: fac(:,:) !!
REAL(r_std), INTENT(in) :: lshead(:,:)
REAL(r_std), INTENT(in) :: hierarchy(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(in) :: fac(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(in) :: lshead(ijdimmax,ijdimmax)
REAL(r_std), DIMENSION(:,:), INTENT(in) :: diaglalo !! Point (in Lat/Lon) where diagnostics will be printed.
!
! Modified
!
INTEGER(i_std), INTENT(inout) :: trip(:,:) !! The trip field (unitless)
INTEGER(i_std), INTENT(inout) :: basin(:,:) !!
REAL(r_std), INTENT(inout) :: topoind(:,:) !! Topographic index of the residence time (m)
INTEGER(i_std), INTENT(inout) :: trip(ijdimmax,ijdimmax) !! The trip field (unitless)
INTEGER(i_std), INTENT(inout) :: basin(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: topoind(ijdimmax,ijdimmax) !! Topographic index of the residence time (m)
!
! For debug and get coordinate of river outlet
!
......@@ -452,6 +453,10 @@ CONTAINS
INTEGER(i_std) :: ie, je !!
INTEGER(i_std) :: sortind(nbvmax) !!
!
INTEGER(i_std) :: rivpas(ijdimmax,ijdimmax) !! Number of passage through grid when following rivers
INTEGER(i_std) :: rivlen(ijdimmax,ijdimmax) !! number of grids along the river
INTEGER(i_std) :: color(ijdimmax,ijdimmax)
!
INTEGER(i_std) :: itrans !!
INTEGER(i_std) :: trans(nbvmax) !!
!
......@@ -975,10 +980,19 @@ CONTAINS
! 7.0 If we calculate topoindex with option 3
! we need to accumulate the topoind_bx along the flow
!
DO ij=1,nb_basin
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
color(ip,jp) = ij
ENDDO
ENDDO
!
! Work only for topo_option number 3 or 4
!
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) THEN
IF ( topo_option .EQ. 3 ) THEN
rivlen(:,:) = 1
rivpas(:,:) = 0
! Loop for all sub-basins and each point belongs to that sub-basin
DO ij=1,nb_basin
DO iz=1,basin_sz(ij)
......@@ -990,7 +1004,8 @@ CONTAINS
p = trip(ip,jp)
cnt = 0
il = ip ; jl = jp
DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj )
! Only follow the river witin the current HTU
DO WHILE ( ij == color(il,jl) .AND. p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj )
cnt = cnt + 1
ill = il + inc(p,1)
jll = jl + inc(p,2)
......@@ -998,18 +1013,38 @@ CONTAINS
! this point still belongs to current sub-basin. Here I
! assume that the sub-routine routing_reg_divbas worked well.
! So we don't need to check.
topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))=topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))+topoind(ill,jll)
topoind(ip,jp)=topoind(ip,jp)+topoind(ill,jll)
rivpas(ill,jll) = rivpas(ill,jll)+1
rivlen(ip,jp) = rivlen(ip,jp)+1
!
il = ill ; jl = jll
p = trip(il,jl)
ENDDO
IF ( topo_option .EQ. 4 ) THEN
! Divide by the number of points of the river to get the mean value.
topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))=topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/cnt
ENDIF
ENDIF
!
ENDDO
! Second pass through all points of the HTU to keep only the topoindex of full rivers.
! These are those throug which we did not pass when following the rivers.
cnt = 0
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
IF ( rivpas(ip,jp) == 0 ) THEN
! Found a river head
topoind(ip,jp) = topoind(ip,jp) * rivlen(ip,jp)
cnt = cnt + rivlen(ip,jp)
ELSE
topoind(ip,jp) = 0.0
ENDIF
ENDDO
! Finish the weight computation of the rivers followed
IF ( cnt > 0 ) THEN
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
topoind(ip,jp) = topoind(ip,jp)/cnt
ENDDO
ENDIF
ENDDO
ENDIF
!
......@@ -1603,9 +1638,9 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!_ ================================================================================================================================
!
!
IF ( topo_option .LT. 1 .OR. topo_option .GT. 4 ) THEN
IF ( topo_option .LT. 1 .OR. topo_option .GT. 3 ) THEN
WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
WRITE(numout,*) ' It should be: 1, 2, 4 or 4'
WRITE(numout,*) ' It should be: 1, 2, or 3'
STOP 'routing_reg_globalize'
ENDIF
!
......@@ -1763,33 +1798,11 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
ENDDO
basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std)
ENDIF
! Minus outlet hierarchy
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4 ) THEN
! All topoindex with HTU have already been weighted and those not at the root of the river set to zero.
IF ( topo_option .EQ. 3 ) THEN
DO iz=1,basin_sz(ij)
area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij))
basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac
basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
ENDDO
!
! Trung: In routing_reg_globalize:
! Trung: Checking negative value of basin_topoind
!
IF ( debug ) THEN
IF ( basin_topoind(ib,ij) .LT. 0. ) THEN
WRITE (numout,*) 'In routing_reg_globalize: We got negative value in ',ib,ij
WRITE (numout,*) 'Start checking ... '
WRITE (numout,*) 'basin_area ', REAL(basin_area(ib,ij))
WRITE (numout,*) ' ==================> '
DO iz=1,basin_sz(ij)
area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij))
WRITE (numout,*) 'iz = ',iz
WRITE (numout,*) 'topoind_bx :', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
WRITE (numout,*) 'area_bx ', area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
WRITE (numout,*) 'area_frac ', area_frac
WRITE (numout,*) 'topoind_bx * area_frac ', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac
ENDDO
ENDIF
ENDIF
!
ENDIF
!
! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy
......@@ -2906,8 +2919,14 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
!basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2)*basin_area(ib, totakeover)&
!& + basin_cg(ib, tokill, 2)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
!!$ IF (outflow_basin(ib,tokill) == totakeover .OR. outflow_basin(ib,totakeover) == tokill ) THEN
!!$ basin_topoind(ib, totakeover) = basin_topoind(ib, totakeover)+basin_topoind(ib, tokill)
!!$ ELSE
!!$ basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
!!$ & + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!!$ ENDIF
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill)
!
......
import sys
import os
import inspect
import importlib
def localdir(follow_symlinks=True):
if getattr(sys, 'frozen', False): # py2exe, PyInstaller, cx_Freeze
path = os.path.abspath(sys.executable)
else:
path = inspect.getabsfile(localdir)
if follow_symlinks:
path = os.path.realpath(path)
return os.path.dirname(path)
def adddirtopath(sdir) :
if not sdir in sys.path :
sys.path.append(localdir()+'/'+sdir)
return
def compilef90(mod, mpi=None) :
adddirtopath(localdir()+"/F90subroutines")
if not "mpi4py" in sys.modules :
err=os.system("cd "+localdir()+"/F90subroutines ; make "+mod)
if err != 0 :
print("Compilation error in the FORTRAN interface")
sys.exit()
else :
mpimod=sys.modules["mpi4py"].MPI
if mpimod.COMM_WORLD.Get_rank() == 0 :
err=os.system("cd "+localdir()+"/F90subroutines ; make "+mod)
if err != 0 :
print("Compilation error in the FORTRAN interface")
sys.exit()
mpimod.COMM_WORLD.Barrier()
return
#
import os
import sys
from inspect import currentframe, getframeinfo
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
#
import numpy as np
from netCDF4 import Dataset
import datetime
......@@ -20,6 +15,10 @@ log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
import F90tools as F90T
#
F90T.adddirtopath("F90subroutines")
F90T.compilef90("routing_interface")
import routing_interface
import NcHelp as NH
import RPPtools as RPP
......@@ -173,7 +172,7 @@ class HydroGraph :
if part.rank == 0:
gvar[gvar >= NCFillValue] = np.nan
#
hist, bin_edges = np.histogram(gvar, bins=nbbins, range=(0,max(10000, np.nanmax(gvar))))
hist, bin_edges = np.histogram(gvar[~np.isnan(gvar)], bins=nbbins, range=(0,min(10000, np.nanmax(gvar))))
topobins = outnf.createVariable("topobins", vtyp, ('nbbnds','nbbins'))
topobins.title = "Topographic index"
topobins.units = "km"
......@@ -184,7 +183,7 @@ class HydroGraph :
topohist.units = "counts"
topohist[:] = hist[:]
#
thist, tbin_edges = np.histogram(gvar*tcst.stream_tcst, bins=nbbins, range=(0,RPP.OneDay/8))
thist, tbin_edges = np.histogram(gvar[~np.isnan(gvar)]*tcst.stream_tcst, bins=nbbins, range=(0,RPP.OneDay/8))
tbins = outnf.createVariable("tstepbins", vtyp, ('nbbnds','nbbins'))
tbins.title = "Time bins"
tbins.units = "s"
......@@ -194,6 +193,12 @@ class HydroGraph :
ncvar.title = "Stable timestep distribution"
ncvar.units = "count"
ncvar[:] = thist[:]
#
maxtstep = np.quantile(gvar[~np.isnan(gvar)], 0.1)
mtstep = outnf.createVariable("MaxTimeStep", vtyp, ('dimpara'), fill_value=NCFillValue)
mtstep.title = "The maximum time-step possible given the distribution of HTU properties"
mtstep.units = "s"
mtstep[:] = maxtstep
return
#
#
......
......@@ -8,20 +8,11 @@ from mpi4py import MPI
import gc
#
import sys
from inspect import currentframe, getframeinfo
#
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
F90=localdir+'/F90subroutines'
if MPI.COMM_WORLD.Get_rank() == 0 :
err=os.system("cd "+F90+"; make all")
if err != 0 :
print("Compilation error in the FORTRAN interface")
sys.exit()
else :
print("Not compiling on other cores")
MPI.COMM_WORLD.Barrier()
import F90tools as F90T
#
F90T.adddirtopath("F90subroutines")
F90T.compilef90("routing_interface")
import routing_interface
import NcHelp as NH
#
......
......@@ -9,17 +9,8 @@ import codecs
import os
from inspect import currentframe, getframeinfo
import sys
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
"""
F90=localdir+'/F90subroutines'
err=os.system("cd "+F90+"; make all")
if err != 0 :
print("Compilation error in the FORTRAN interface")
sys.exit()
"""
from diagst import *
#
import F90tools as F90T
#
import getargs
config = getargs.SetupConfig()
......@@ -28,6 +19,10 @@ log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
F90T.adddirtopath("F90subroutines")
F90T.compilef90("diagst")
import diagst
#
class Locations :
def __init__(self, bbox) :
import RPPtools as RPP
......@@ -279,6 +274,7 @@ class Locations :
#
class PlacedLocations :
def __init__(self, GraphFile) :
#
nc = Dataset(GraphFile, 'r')
self.nbloc = nc.dimensions["locations"].size
self.gridshape = [nc.dimensions["y"].size,nc.dimensions["x"].size]
......@@ -331,7 +327,7 @@ class PlacedLocations :
"""
g0 = self.nbpt_glo[int(j)-1,int(i)-1] # -1 Because F -> C
b0 = int(k)
mask = diagst.upstmask(
mask = diagst.diagst.upstmask(
nbpt = self.nbpt,
nbasmax = self.nbasmax,
routetogrid = self.routetogrid,
......
......@@ -6,10 +6,6 @@ import os, sys
from mpi4py import MPI
import gc
#
import matplotlib as mpl
mpl.use('Agg')
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import time
#
# Version of RoutingPP
......
......@@ -7,20 +7,8 @@ from inspect import currentframe, getframeinfo, getabsfile
from netCDF4 import Dataset
from netCDF4 import stringtochar, chartostring, stringtoarr
#
def localdir(follow_symlinks=True):
if getattr(sys, 'frozen', False): # py2exe, PyInstaller, cx_Freeze
path = os.path.abspath(sys.executable)
else:
path = getabsfile(localdir)
if follow_symlinks:
path = os.path.realpath(path)
return os.path.dirname(path)
localdir=localdir()
#localdir=os.path.dirname(getframeinfo(currentframe()).filename)
#localdir=os.path.dirname(os.path.realpath(__file__))
sys.path.append(localdir+'/Stations')
sys.path.append(localdir+'/F90subroutines')
import F90tools as F90T
F90T.adddirtopath("Stations")
FillValue = 1.e+20
#
# Gather user input.
......@@ -66,20 +54,7 @@ A.fetchnetcdf(GRDCFile, metaonly=True, region=[L.lonrange,L.latrange])
#
IDs=[]
initial_calls=[]
#create a FindID function that calls other functions depending on the input
for teststn in StationList :
#next step: when there is a @ in the string look at what is before
#if there is a name!=all it is a river, if it's a nuber it's a region or subregion
#this info must go in input to the ID finder
#next next step: add an exclusion list that uses the same function as below but creates a list of ID_excluded to subtract from the IDs list
#for the moment just add:
#elif there is a @ in the string:
#if after @ there is an integer
#if <10 -> region -> FindIDfromWMOreg and extend
#else -> subregion -> FindIDfromWMOsubreg and extend
#elif after @ is a string
#if 2 char and isupper true -> country -> FindIDfromCountry and extend
#else -> river -> FindIDfromRiver and extend
IDs.extend(A.FindID(teststn))
for i in A.FindID(teststn):
initial_calls.append(teststn)
......
......@@ -11,19 +11,20 @@ import numpy as np
import numpy.ma as ma
import time
import sys
import os
from inspect import currentframe, getframeinfo
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
import routing_interface
#
import getargs
config = getargs.SetupConfig()
frac_lim=config.getfloat("OverAll", "frac_lim", fallback=0.3)
limit_step5=config.getfloat("OverAll","limit_step5",fallback=5)
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
#
import F90tools as F90T
#
F90T.adddirtopath("F90subroutines")
F90T.compilef90("routing_interface")
import routing_interface
#
######################
class trunc:
......
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
import getargs
log_master, log_world = getargs.getLogger()
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
inc=np.zeros((8,2), dtype=np.int8)
inc[0,0] = 1; inc[0,1] =0
inc[1,0] = 1; inc[1,1] =1
inc[2,0] = 0; inc[2,1] =1
inc[3,0] = -1; inc[3,1] =1
inc[4,0] = -1; inc[4,1] =0
inc[5,0] = -1; inc[5,1] =-1
inc[6,0] = 0; inc[6,1] =-1
inc[7,0] = 1; inc[7,1] =-1
epsilon=0.00001
undef_int=999999999
#
def buildfullgrid(lon, lat, lonint, latint) :
nbpts, nlon, nlat = lon.shape
# LON order
for i in range(nlat) :
if np.count_nonzero(~np.isnan(lon[0,:,i])) > 0 :
dlon = np.sign(np.gradient(lon[0,:,i]))[~np.isnan(np.gradient(lon[0,:,i]))][0]
# LAT order
for i in range(nlon) :
if np.count_nonzero(~np.isnan(lat[0,i,:])) > 0 :
dlat = np.sign(np.gradient(lat[0,i,:]))[~np.isnan(np.gradient(lat[0,i,:]))][0]
#
xlon=np.unique(np.hstack(lon))[~np.isnan(np.unique(np.hstack(lon)))]
xlat=np.unique(np.hstack(lat))[~np.isnan(np.unique(np.hstack(lat)))]
# Zomm in
xlon=xlon[np.argmin(np.abs(xlon-min(lonint))):np.argmin(np.abs(xlon-max(lonint)))]
xlat=xlat[np.argmin(np.abs(xlat-min(latint))):np.argmin(np.abs(xlat-max(latint)))]
#
if np.abs(np.min(np.gradient(xlon))-np.min(np.gradient(xlon))) < epsilon :
res_lon = np.mean(np.gradient(xlon))
else :
ERROR("Could not determine the longitude resolution of the hydrological grid")
if np.abs(np.min(np.gradient(xlat))-np.min(np.gradient(xlat))) < epsilon :
res_lat = np.mean(np.gradient(xlat))
else :
ERROR("Could not determine the latitude resolution of the hydrological grid")
if np.sign(np.gradient(xlon))[0] != dlon :
lonv=np.concatenate(([np.max(xlon)+res_lon],xlon[::-1],[np.min(xlon)-res_lon]))
else :
lonv=np.concatenate(([np.min(xlon)-res_lon],xlon[:],[np.max(xlon)+res_lon]))
if np.sign(np.gradient(xlat))[0] != dlat :
latv=np.concatenate(([np.max(xlat)+res_lat],xlat[::-1],[np.min(xlat)-res_lat]))
else :
latv=np.concatenate(([np.min(xlat)-res_lat],xlat[:],[np.max(xlat)+res_lat]))
#
# return the full grid where longitude and latitude are 2D arrays
#
return np.meshgrid(lonv, latv, indexing='ij')
#
#
#
def enlargegrid(lon, lat, trip, data, lon_full, lat_full) :
#
nbpts, nlon, nlat = lon.shape
nlonfull,nlatfull= lon_full.shape
fulldata=np.zeros((nlonfull,nlatfull), dtype=np.float32)
fulldata[:,:]=np.nan
fulltrip=np.zeros((nlonfull,nlatfull), dtype=np.float32)
fulltrip[:,:]=undef_int
for ib in range(nbpts) :
for i in range(nlon) :
for j in range(nlat) :
if trip[ib,i,j] < 1000 :
dist=np.sqrt((lon_full-lon[ib,i,j])**2 + (lat_full-lat[ib,i,j])**2)
il,jl=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
fulltrip[il,jl]=trip[ib,i,j]
fulldata[il,jl]=data[ib,i,j]
return fulltrip, fulldata
#
#
#
def enlargebasins(lon, lat, datasz, data, lon_full, lat_full) :
nbpts, nlon, nlat = lon.shape
nlonfull, nlatfull = lon_full.shape
nbpts, nbvmax = datasz.shape
fullbasins=np.zeros((nlonfull,nlatfull), dtype=np.float32)
fullbasins[:,:]=np.nan
for ib in range(nbpts) :
for iz in range(nbvmax) :
if datasz[ib,iz] > 0 :
for isz in range(datasz[ib,iz]) :
# Fortran to C indexing
i=data[ib,iz,isz,0]-1
j=data[ib,iz,isz,1]-1
dist=np.sqrt((lon_full-lon[ib,i,j])**2 + (lat_full-lat[ib,i,j])**2)
il,jl=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
fullbasins[il,jl]=iz
return fullbasins
#
#
#
def closestpoint(lonout, latout, nbptb, basin_pts, lonbx, latbx) :
lonin = 3.5
latin = 39.0
dist = undef_int