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

Commit be3594b0 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Adds the bounds to the graph file for all grid points and not only the land...

Adds the bounds to the graph file for all grid points and not only the land points. This will facilitate the interpolations with CDO.
parent 4e8c98c1
......@@ -69,6 +69,7 @@ def addparameters(outnf, version, part, tcst, vtyp, NCFillValue, lim_floodcri) :
# Add global attributes to the netCDF file
outnf.setncattr("RoutingPPVersion", version)
outnf.setncattr("GenerationDate", datetime.datetime.utcnow().isoformat(sep=" "))
outnf.setncattr("Conventions", "CF-1.6")
# Add routing parameters
outnf.createDimension('dimpara', 1)
stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
......@@ -293,7 +294,6 @@ class HydroGraph :
#
vtyp=np.float64
cornerind=[0,2,4,6]
nbcorners = len(cornerind)
nbbins = 2000
#
nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
......@@ -305,7 +305,7 @@ class HydroGraph :
outnf.createDimension('y', globalgrid.nj)
outnf.createDimension('land', len(procgrid.area))
outnf.createDimension('z', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('bnd', len(cornerind))
outnf.createDimension('inflow', self.max_inflow)
outnf.createDimension('stnperhtu', nlmax)
outnf.createDimension('nbbins', nbbins)
......@@ -314,7 +314,7 @@ class HydroGraph :
else :
outnf = None
#
NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue)
NH.addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
addparameters(outnf, version, part, tcst, vtyp, NCFillValue, self.lim_floodcri)
#
......
......@@ -31,16 +31,22 @@ def mindist(coord,lon,lat) :
d.append(RPP.loladist(np.array(c), np.array([lon,lat])))
return np.argmin(d)
#
#
#
def genindex(ni, nj, indFi, indFj, land=None) :
if land is None :
indP = [[j,i] for i in range(ni) for j in range(nj)]
else :
indP = [[j,i] for i in range(ni) for j in range(nj) if land[j,i] > 0]
indF = [[indFi[j,i],indFj[j,i]] for j,i in indP]
return indP, indF
#
# Function to gather all land points but while keeping also the neighbour information.
#
def gatherland(lon, lat, land, indP, indFi, indFj) :
def gatherland(indP, lon, lat, land) :
nj,ni=lon.shape
#
indP_land = [[j,i] for i in range(ni) for j in range(nj) if land[j,i] > 0]
coord = [[lon[j,i],lat[j,i]] for j,i in indP_land]
indF_land = [[indFi[j,i],indFj[j,i]] for j,i in indP_land]
#
nbland=len(coord)
coord = [[lon[j,i],lat[j,i]] for j,i in indP]
#
# Get the neighbours in the coord list. The same order is used as in ORCHIDEE
#
......@@ -49,9 +55,9 @@ def gatherland(lon, lat, land, indP, indFi, indFj) :
# -1 will become 0 and indicate point is outside of domain.
# -2 will become -1 and indicate ocean.
#
neighbours = [get_neighbours(i,j,ni,nj,land, coord, lon, lat) for j,i in indP_land]
neighbours = [get_neighbours(i, j, ni, nj, land, coord, lon, lat) for j,i in indP]
#
return nbland, coord, neighbours, indP_land, indF_land
return coord, neighbours
#
def get_neighbours(i,j, ni, nj, land, coord, lon, lat):
nn = [[j+r[0], i +r[1]] for r in rose]
......@@ -59,8 +65,9 @@ def get_neighbours(i,j, ni, nj, land, coord, lon, lat):
return ntmp
#
#
# Generate lists (corners and center) of Lat/Lon polygones for the mesh.
#
def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
def lalopoly(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
if lon.shape[1]>1 and lat.shape[0]>1:
dlon=int((lon[0,-1]-lon[0,0])/np.abs(lon[0,-1]-lon[0,0]))
......@@ -74,6 +81,12 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
Lpolyg = [RPP.boxit(np.array([istart+ij[0],jstart+ij[1]]), dlon, dlat, 2) for ij in indF]
cornersll = [proj.ijll(polyg) for polyg in Lpolyg]
centersll = [proj.ijll([[istart+ij[0],jstart+ij[1]]])[0] for ij in indF]
#
return cornersll, centersll
#
# Transfer the Lat/Lon polygones to SphericaGeometry polygones
#
def sphericalpoly(cornersll, centersll) :
#
# Avoid that points have longitude or latitude exactly equal to zero
#
......@@ -87,7 +100,7 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
box=[[np.min(np.array(allon)),np.max(np.array(allon))],[np.min(np.array(allat)),np.max(np.array(allat))]]
#
return cornersll, cornerspoly, centersll, radiusll, areas, box
return cornerspoly, radiusll, areas, box
#
# Extract the coordinates of the grid file. If ni < 0 and nj < 0 then the
# full grid is read.
......@@ -262,11 +275,20 @@ class ModelGrid :
indFi=[[i+1] for j in range(nj) for i in range(ni)]
indFj=[[j+1] for j in range(nj) for i in range(ni)]
#
# Gather all land points
# Build index lists for all points and land points
#
self.indP_full, self.indF_full = genindex(ni, nj, np.reshape(indFi[:],self.lon_full.shape),\
np.reshape(indFj[:],self.lon_full.shape))
self.indP, self.indF = genindex(ni, nj, np.reshape(indFi[:],self.lon_full.shape),\
np.reshape(indFj[:],self.lon_full.shape), land=self.land)
#
# Build lists of central coordinates and neighbours of land points
#
self.coordll,self.neighbours = gatherland(self.indP, self.lon_full, self.lat_full, self.land)
self.nbland = len(self.coordll)
#
#
#
self.nbland, self.coordll,self.neighbours,self.indP,indF = gatherland(self.lon_full,self.lat_full,self.land,ind,\
np.reshape(indFi[:],self.lon_full.shape),\
np.reshape(indFj[:],self.lon_full.shape))
INFO("Shape of region :"+str(ni)+" x "+str(nj)+" with nbland="+str(self.nbland))
if self.nbland < 1 or self.nbland > self.nj*self.ni :
INFO("Found impossible number of land points : "+str(self.nbland))
......@@ -279,20 +301,20 @@ class ModelGrid :
#
for ip in self.neighbours[0][:] :
if ip >= 0 :
DEBUG("Neighbour : "+str(ip)+" P index : "+str(self.indP[ip])+" F index : "+str(indF[ip][:]))
DEBUG("Neighbour : "+str(ip)+" P index : "+str(self.indP[ip])+" F index : "+str(self.indF[ip][:]))
else :
DEBUG("Neighbour : "+str(ip)+" Not Land")
#
# Define projection
#
if griddesc['type'] == "lambert_conformal_conic" :
proj=PS.LambertC(griddesc['dx'], griddesc['known_lon'], griddesc['known_lat'], griddesc['truelat1'], \
self.proj=PS.LambertC(griddesc['dx'], griddesc['known_lon'], griddesc['known_lat'], griddesc['truelat1'], \
griddesc['truelat2'], griddesc['stdlon'])
elif griddesc['type'] == "rotated_pole" :
proj = PS.Cassini(griddesc['lon0'], griddesc['lat0'], griddesc['stdlon'], griddesc['lon1'], griddesc['lat1'], \
self.proj = PS.Cassini(griddesc['lon0'], griddesc['lat0'], griddesc['stdlon'], griddesc['lon1'], griddesc['lat1'], \
griddesc['dlon'], griddesc['dlat'])
elif griddesc['type'] == "RegLonLat" :
proj=PS.RegLonLat(self.res_lon, self.res_lat, griddesc['inilon'], griddesc['inilat'])
self.proj=PS.RegLonLat(self.res_lon, self.res_lat, griddesc['inilon'], griddesc['inilat'])
else :
ERROR("Unknown grid type")
sys.exit()
......@@ -301,8 +323,8 @@ class ModelGrid :
#
self.box=[[np.min(self.lon_full),np.max(self.lon_full)],[np.min(self.lat_full),np.max(self.lat_full)]]
#
self.polyll, self.polylist, self.centers, self.radius, self.area, self.box_land = \
corners(indF, proj, istart, jstart, self.lon_full, self.lat_full, self.res_lon, self.res_lat)
self.polyll, self.centers = lalopoly(self.indF, self.proj, istart, jstart, self.lon_full, self.lat_full, self.res_lon, self.res_lat)
self.polylist, self.radius, self.area, self.box_land = sphericalpoly(self.polyll, self.centers)
#
geo.close()
#
......@@ -451,9 +473,27 @@ class GlobalGrid :
self.latmat = lat_full[self.jgstart:self.jgstart+self.nj,self.igstart:self.igstart+self.ni]
INFO("Shape of full grid kept in GlobalGrid : "+str(self.lonmat.shape))
if self.lonmat.shape[1]>1 and self.latmat.shape[0]>1:
self.londir=int((self.lonmat[0,-1]-self.lonmat[0,0])/np.abs(self.lonmat[0,-1]-self.lonmat[0,0]))
self.latdir=int((self.latmat[0,0]-self.latmat[-1,0])/np.abs(self.latmat[0,0]-self.latmat[-1,0]))
else :
self.londir = res_lon
self.latdir = res_lat
self.land = getland(geo, self.igstart, self.ni, self.jgstart, self.nj)
self.nbland = int(np.sum(self.land))
#
geo.close()
#
# Function to add the corners to the grid
#
def updatecorners(self, modelgrid) :
cornerind=[0,2,4,6]
self.lonmat_bnd = np.zeros((len(cornerind),self.nj,self.ni))
self.latmat_bnd = np.zeros((len(cornerind),self.nj,self.ni))
for j in range(self.nj) :
for i in range(self.ni) :
corners = modelgrid.proj.ijll(RPP.boxit(np.array([self.igstart+i+1,self.jgstart+j+1]), self.londir, self.latdir, 2))
self.lonmat_bnd[:,j,i] = [corners[ic][0] for ic in cornerind]
self.latmat_bnd[:,j,i] = [corners[ic][1] for ic in cornerind]
return
......@@ -2,65 +2,64 @@
import numpy as np
#
#
def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind) :
def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue) :
#
# Longitude
longitude = part.gather(procgrid.lon_full)
if part.rank == 0 :
lon=outnf.createVariable("lon", vtyp, ('y','x'), fill_value=NCFillValue)
lon.units="grid box centre degrees east"
lon.title="Longitude"
lon.units="degrees_east"
lon.long_name="Longitude"
lon.standard_name="grid_longitude"
lon.axis="X"
lon._CoordinateAxisType = "Lon"
lon.bounds="lon_bnd"
lon[:,:] = globalgrid.lonmat[:,:]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
if part.rank == 0 :
lat=outnf.createVariable("lat", vtyp, ('y','x'), fill_value=NCFillValue)
lat.units="grid box centre degrees north"
lat.standard_name="grid latitude"
lat.title="Latitude"
lat.units="degrees_north"
lat.standard_name="grid_latitude"
lat.long_name="Latitude"
lat.axis="Y"
lat._CoordinateAxisType = "Lat"
lat.bounds="lat_bnd"
lat[:,:] = globalgrid.latmat[:,:]
#
# Bounds of grid box
#
llonpoly=np.zeros((nbcorners,procgrid.nbland))
llatpoly=np.zeros((nbcorners,procgrid.nbland))
for i in range(procgrid.nbland) :
llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
lon_bnd = procgrid.landscatter(llonpoly)
lat_bnd = procgrid.landscatter(llatpoly)
if part.rank == 0 :
lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
lonbnd.units="grid box corners degrees east"
lonbnd.units="degrees_east"
lonbnd.title="Longitude of Corners"
latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
latbnd.units="grid box corners degrees north"
latbnd.units="degrees_north"
latbnd.title="Latitude of Corners"
lonbnd[:,:,:] = globalgrid.lonmat_bnd[:,:,:]
latbnd[:,:,:] = globalgrid.latmat_bnd[:,:,:]
else :
lonbnd= np.zeros((1,1,1))
latbnd= np.zeros((1,1,1))
lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
#
# Land sea mask
#
if part.rank == 0 :
land=outnf.createVariable("land", vtyp, ('y','x'), fill_value=NCFillValue)
land.units="Land Sea mask"
land.standard_name="landsea mask"
land.long_name="landsea mask"
land.standard_name="land_area_fraction"
land.title="Land"
land.coordinates="lat lon"
land[:,:] = globalgrid.land[:,:]
# Area
areas = procgrid.landscatter(np.array(procgrid.area, dtype=np.float64))
areas[np.isnan(areas)] = NCFillValue
if part.rank == 0 :
area=outnf.createVariable("area", vtyp, ('y','x'), fill_value=NCFillValue)
area=outnf.createVariable("Areas", vtyp, ('y','x'), fill_value=NCFillValue)
area.units="m^2"
area.standard_name="grid area"
area.long_name="Mesh areas"
area.title="Area"
area.coordinates="lat lon"
else :
area = np.zeros((1,1))
area[:,:] = part.gather(areas[:,:])
......
......@@ -53,6 +53,7 @@ rank=comm.Get_rank()
part = PA.partition(gg.ni, gg.nj, gg.land, comm, nbcore, szhalo, rank)
#
modelgrid=MG.ModelGrid(part.ihstart+gg.igstart,part.nih,part.jhstart+gg.jgstart,part.njh)
gg.updatecorners(modelgrid)
#
INFO("Longitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[0][0])+" : "+str(modelgrid.box_land[0][1]))
INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+" : "+str(modelgrid.box_land[1][1]))
......
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