Commit 6a89d8a3 authored by Anthony's avatar Anthony
Browse files

Tools : HydroDATA_shapefile

parent 2ef3ca4a
......@@ -39,7 +39,7 @@ xcuserdata/
Weights/
tests/*/Weights/
tests/*/DocumentationInterface
Tools/HydroDATA_shapefile/shapefile/*
Tools/HydroDATA_shapefile/shapefile/*/
Tools/HydroDATA_shapefile/*.nc
__pycache__/
./run.def
......
#!/bin/bash
#
#########################################
#
# Script to build and maintain en environment to run
# the python code in this directory. We verify that
# all needed packages are available.
# This is based on the anaconda python system : https://www.anaconda.com/distribution/
#
# Usage : source Environment
#
#########################################
#
# Functions
#
function hasanaconda (){
ipsl="ipsl camelot climserv ciclad"
idris="idris jean-zay"
case $(hostname -d) in
climserv*|private.ipsl.fr|ipsl.polytechnique.fr)
echo "IPSL"
eval "$1='ipsl'"
;;
idris*|jean-zay*)
echo "IDRIS"
eval "$1='idris'"
;;
*)
echo "unknown machine : $(hostname -d)"
eval "$1='unknown'"
esac
}
function checkenv (){
rm -f checkenv.py
cat <<EOF > checkenv.py
import mpi4py
import spherical_geometry
import netCDF4
import matplotlib
import mpl_toolkits
import configparser
import numba
import cartopy
print ("OK")
EOF
python checkenv.py > /dev/null
if [ $? -gt 0 ] ; then
echo "Environment not OK. Some packages must be missing."
echo "This was tested with the checkenv.py code."
else
rm -f checkenv.py
fi
}
#
# Main
#
loc=''
hasanaconda loc
if [ loc != "unknown" ] ; then
#
# Make sure the right version of Python
#
case $loc in
"ipsl")
module unload python
module load python/3.6-anaconda50
#
# Test if the MPI environment is installed.
#
if [ ! -e ${HOME}/.conda/envs/MPI ] ; then
#
# Create the MPI environment as it does not exist.
#
conda create -y -n MPI mpi4py numba astropy netcdf4 matplotlib basemap ConfigParser
source activate MPI
conda install -y cartopy
conda install -y --channel astropy spherical-geometry
else
source activate MPI
fi
;;
"idris")
module load python/3.7.3
if [[ -L ${HOME}/.conda && -d ${HOME}/.conda ]] ; then
if [ ! -e ${WORK}/.conda/envs/MPI ] ; then
conda create -y -n MPI mpi4py numba astropy netcdf4 matplotlib basemap ConfigParser cartopy
conda activate MPI
if [[ -L ${HOME}/.local && -d ${HOME}/.local ]] ; then
export PYTHONUSERBASE=$WORK/.local
if [ $(find $PYTHONUSERBASE -name spherical_geometry -ls | wc -l) -le 0 ] ; then
git clone https://github.com/spacetelescope/spherical_geometry.git
cd spherical_geometry
python setup.py install --user
fi
else
echo "Make sure \$HOME/.link is a link to your \$WORK/.local"
fi
else
conda activate MPI
fi
else
echo "Make sure \$HOME/.conda is a link to your \$WORK/.conda"
fi
;;
*)
echo "Environment for $loc is not foreseen yet"
exit
;;
esac
fi
#
# Verify we can load all the needed modules.
#
checkenv
#!/bin/bash
#
#PBS -N remap_MERIT
#
#PBS -j oe
#PBS -l nodes=2:ppn=32
#PBS -l walltime=128:00:00
#PBS -l vmem=300gb,mem=300gb
#
cd ${PBS_O_WORKDIR}
#
# Set the right Python 3 Anaconda environment
#
source ./Environment
#
\rm *.txt
python ./main_multiprocess.py
ls -l
from netCDF4 import Dataset
import numpy as np
from numba import jit
#from spherical_geometry import polygon
from shapely.geometry import Polygon
import time
import sys
import ModelGrid as MG
import RPPtools as RPP
import Projections as PS
import cartopy.io.shapereader as shpreader
from shapely.geometry import MultiPolygon, Polygon
EarthRadius=6370000.0
######################################
class partition:
def __init__(self, nbcore, nloop, numpolyg):
self.nbcore = nbcore
self.nbloop = nloop
self.numpolyg = numpolyg
n = int(numpolyg / (nbcore*nloop))
self.inter = [[n*i,n*(i+1)] for i in range(nbcore*nloop-1)]
self.inter.append([n*(nbcore*nloop-1),numpolyg])
def get_part(self, i):
return self.inter[i]
#####################################
class GlobalGrid:
def __init__(self, fdir, input_type):
# Get the coordinates
nc = Dataset(fdir, "r")
griddesc, self.lon, self.lat, self.res_lon, self.res_lat = MG.getcoordinates(nc, 0, -1, 0, -1)
self.nj,self.ni = self.lon.shape
#
self.land = np.full(self.lon.shape, 1); self.nbland = np.sum(self.land)
if input_type == "geogrid":
# Get the projection
if griddesc['type'] == "lambert_conformal_conic" :
self.proj = PS.LambertC(griddesc['dx'], griddesc['known_lon'], griddesc['known_lat'], griddesc['truelat1'], \
griddesc['truelat2'], griddesc['stdlon'])
elif griddesc['type'] == "rotated_pole" :
self.proj = PS.Cassini(griddesc['lon0'], griddesc['lat0'], griddesc['stdlon'], griddesc['lon1'], griddesc['lat1'], \
griddesc['dlon'], griddesc['dlat'])
elif griddesc['type'] == "RegLonLat" :
self.proj=PS.RegLonLat(self.res_lon, self.res_lat, griddesc['inilon'], griddesc['inilat'])
else :
ERROR("Unknown grid type")
sys.exit()
elif input_grid == "hydro":
try:
self.proj=PS.RegLonLat(self.res_lon, self.res_lat, self.lon[0], self.lat[0])
except:
print("Error in Hydrogrid projection")
print("res_lon:{0}, res_lat:{1}, lon0: {2}, lat0: {3}".format(self.res_lon, self.res_lat, self.lon[0], self.lat[0]))
sys.exit()
self.grid_params()
self.get_polygons()
def grid_params(self):
self.igstart = 0
self.jgstart = 0
# Change name self.ind
indFi=[[i+1] for j in range(self.nj) for i in range(self.ni)]
indFj=[[j+1] for j in range(self.nj) for i in range(self.ni)]
indP_land = [[j,i] for i in range(self.ni) for j in range(self.nj)]
indFi=np.reshape(indFi, self.lon.shape)
indFj=np.reshape(indFj, self.lon.shape)
self.indF = [[indFi[j,i],indFj[j,i]] for j,i in indP_land]
#self.coord = [[self.lon[j,i],self.lat[j,i]] for i,j in self.indF]
#
self.dlon=int((self.lon[0,-1]-self.lon[0,0])/np.abs(self.lon[0,-1]-self.lon[0,0]))
self.dlat=int((self.lat[0,0]-self.lat[-1,0])/np.abs(self.lat[0,0]-self.lat[-1,0]))
def get_polygons(self):
print("BORDERS")
t0 = time.time()
Lpolyg = [RPP.boxit(np.array([self.igstart + ij[0],self.jgstart + ij[1] ]), self.dlon, self.dlat, 2) for ij in self.indF]
t1 = time.time()
print("{0} s.".format(int(t1-t0)))
#
print("Corners")
t0 = time.time()
cornersll = [self.proj.ijll(polyg) for polyg in Lpolyg]
t1 = time.time()
print("{0} s.".format(int(t1-t0)))
#
print("Poly")
t0 = time.time()
self.polygg = [Polygon(c) for c in cornersll]
t1 = time.time()
print("{0} s.".format(int(t1-t0)))
######################################
class landPolygon:
def __init__(self, land_dir):
print("Load Land")
t0 = time.time()
shpfi = shpreader.Reader(land_dir)
for k, r in enumerate(shpfi.geometries()):
if k == 0:
self.P = MultiPolygon(r)
else:
self.P.union(MultiPolygon(r))
t1 = time.time()
print("{0} s.".format(int(t1-t0)))
import numpy as np
import numpy.ma as ma
import time
#
EarthRadius=6370000.0
#
def interpolate(M, F, option = "conservative"):
anc_new = ma.zeros(M.trip.shape)
anc_new.mask = M.trip.mask
for c, value in zip(F.centers, F.values):
# Polygone from floodplains
hydrocell = F.corners(c)
r = F.hd
# Corresponding MERIT grid points
indices, polylist, subcenters = M.select(c,r)
# F.fp[jfp, ifp] = 1 si sélectionné
OVERLAP = []
if len(indices) > 0:
for i, cell in enumerate(polylist):
try:
area_in = cell.intersection(hydrocell).area
except:
area_in = 0
if area_in > 0:
frac_modelgrid = area_in/cell.area
frac_input = area_in/hydrocell.area
OVERLAP.append([indices[i], frac_modelgrid, frac_input])
# Now distribute the value in function of the chosen option
# The overlap may be the opportunity ti distribute with
# an order of priority
# Ex : value of input is less than 100% and we want to distribute it
# over the lowest model grid points
if option == "conservative":
for ind, frac_modelgrid, frac_input in OVERLAP:
jm, im = M.ind_land[ind]
anc_new[jm,im] = anc_new[jm,im] + frac_modelgrid*value
return anc_new
#
from netCDF4 import Dataset
import numpy as np
import numpy.ma as ma
def create_output(dinput, P_list,P, gg, output_name):
#
nc = Dataset(dinput, "r")
lon = nc.variables["XLONG_M"][0,:,:]
lat = nc.variables["XLAT_M"][0,:,:]
njg, nig = gg.lon.shape
# Get the data
outland = np.zeros((njg, nig))
# Mask the output
# it will be unmasked where we processed the ancillary data
#
# P is partition
for k, out_part in enumerate(P_list):
i0, i1 = P.get_part(k)
for l, m in zip(range(i0,i1),range(len(out_part))):
i,j = gg.indF[l]
outland[j-1,i-1] = out_part[m]
#ancnew = ma.masked_where(ancnew == -1, ancnew)
# Start the netCDF
with Dataset(output_name,'w') as foo:
dlat = foo.createDimension('y', (njg))
dlon = foo.createDimension('x', (nig))
#
latitudes = foo.createVariable('lat', np.float32,('y','x'))
longitudes = foo.createVariable('lon', np.float32,('y','x'))
#
latitudes[:,:] = gg.lat[:,:]
longitudes[:,:] = gg.lon[:,:]
#
# Fonction de name_var !!
out = foo.createVariable("landfrac", np.float32,('y','x'))#, zlib = True)
out[:,:] = outland[:,:]
foo.sync()
import sys
sys.path.append("./lib")
sys.path.append("../../")
import numpy as np
from netCDF4 import Dataset
import multiprocessing
#############
ConfigFile="run.def"
import getargs
config = getargs.SetupConfig()
Debug = config.getboolean("OverAll", "Debug", fallback = False)
dinput=config.get("OverAll", "dinput", fallback = None)
shpdir=config.get("OverAll", "shpdir", fallback = None)
input_type=config.get("OverAll", "input_type", fallback = None)
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
#
nloop=config.getint("OverAll", "nloop", fallback=100)
nbcore=config.getint("OverAll", "nbcore", fallback=1)
#
doutput=config.get("OverAll", "doutput", fallback = None)
#
from GRID import GlobalGrid, landPolygon, partition
from output import create_output
#############
if (input_type is None):
print("Please define input_type")
sys.exit()
elif (input_type not in ["geogrid", "hydro"]):
print("input_type= {0} is not a valid choiche -> = geogrid or = hydro".format(input_type))
elif (dinput is None):
print("Please define dinput")
sys.exit()
elif (shpdir is None):
print("Please define shpdir")
sys.exit()
elif (doutput is None):
print("Please define doutput")
sys.exit()
if multiprocessing.cpu_count() < nbcore:
print("You only have {0} CPUs available".format(multiprocessing.cpu_count()))
print("Please reduce nbcore (={0}) or raise the number of CPUs available".format(nbcore))
##############
gg = GlobalGrid(dinput, "geogrid")
landP= landPolygon(shpdir)
#############
# MAIN FUNCTION FOR PARALLELIZATION
def calculate(part_rank, test = Debug):
i0,i1 = P.get_part(part_rank)
L = []
for k in range(i0,i1):
if test :
L.append(1)
else:
L.append(gg.polygg[k].intersection(landP.P).area/gg.polygg[k].area)
with open("./out.txt", "a") as foo:
foo.write("part_rank: {0} \n".format(part_rank))
return L
############
# Prepare the partition
P = partition(nbcore,nloop,len(gg.indF))
# Prepare parallelization pool
pool = multiprocessing.Pool(processes=nbcore)
# Launch parallelization
P_list = pool.map(calculate, np.arange(nbcore* nloop)) # nbcore * nloop
############
# Create the output dataset
create_output(dinput, P_list, P,gg, doutput)
[OverAll]
#
EarthRadius = 6370000.
#
# nbcore : number of CPUs to be used
# nbloo : number of loop -> Decomposition in (nbcore * nloop) substeps
# These steps are solved through competitive parralelism
#
nloop = 100
nbcore = 64
# If Debug = True, there will not have overlap calculation (everything = 1)
# Quicker because it s the heavier step, allow to check if everything works
#
Debug = False
# [DEF] input_type : geogrid, hydro
#
input_type =
dinput =
# output file direction
doutput = ./output.nc
# Shapefile direction
shpdir =
[OverAll]
#
EarthRadius = 6370000.
#
# nbcore : number of CPUs to be used
# nbloo : number of loop -> Decomposition in (nbcore * nloop) substeps
# These steps are solved through competitive parralelism
#
nloop = 100
nbcore = 64
# If Debug = True, there will not have overlap calculation (everything = 1)
# Quicker because it s the heavier step, allow to check if everything works
#
Debug = False
# [DEF] input_type : geogrid, hydro
#
input_type = geogrid
dinput = /bdd/MEDI/workspaces/polcher/WRF_Forcing/AmSud-A/geo_em.d01.nc
# output file direction
doutput = ./output.nc
# Shapefile direction
shpdir = ./shapefile/ne_10m_land/ne_10m_land.shp
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment