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

Adding the rotated pole (Cassini) projection to the code.

parent b6d25b64
......@@ -22,6 +22,7 @@ INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
rose=[[-1,-1],[-1,0],[-1,+1],[0,+1],[+1,+1],[+1,0],[+1,-1],[0,-1]]
epsilon=0.00001
projcodes={"rotated_pole":6,"latitude_longitude":5,"lambert_conformal_conic":1}
#
###################################################################################
#
......@@ -99,8 +100,11 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
# Guess the type of grid
#
griddesc = {}
if "DX" in geo.ncattrs() :
griddesc['type'] = "RegXY"
if "MAP_PROJ" in geo.ncattrs() :
code=geo.MAP_PROJ
for p in projcodes.keys() :
if projcodes[p] == code :
griddesc['type'] = p
elif len(geo.variables["lon"].shape) == 1 :
griddesc['type'] = "RegLonLat"
else :
......@@ -109,18 +113,39 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
#
# Extract grid information
#
if griddesc['type'] == "RegXY" :
# We have a geogrid file from WRF
griddesc['dx'] = geo.DX
griddesc['known_lon'] = geo.corner_lons[0]
griddesc['known_lat'] = geo.corner_lats[0]
griddesc['truelat1'] = geo.TRUELAT1
griddesc['truelat2'] = geo.TRUELAT2
griddesc['stdlon'] = geo.STAND_LON
if griddesc['type'] != "RegLonLat" :
# We have a geogrid file from which to get information
if griddesc['type'] == "lambert_conformal_conic" :
griddesc['dx'] = geo.DX
griddesc['known_lon'] = geo.corner_lons[0]
griddesc['known_lat'] = geo.corner_lats[0]
griddesc['truelat1'] = geo.TRUELAT1
griddesc['truelat2'] = geo.TRUELAT2
griddesc['stdlon'] = geo.STAND_LON
elif griddesc['type'] == "rotated_pole" :
griddesc['stdlon'] = geo.STAND_LON
griddesc['lon0'] = geo.POLE_LON
griddesc['lat0'] = geo.POLE_LAT
griddesc['known_lon'] = geo.corner_lons[0]
griddesc['known_lat'] = geo.corner_lats[0]
griddesc['dlon'] = geo.DX/(EarthRadius * np.pi * 2.0 / 360.0)
griddesc['dlat'] = geo.DY/(EarthRadius * np.pi * 2.0 / 360.0)
griddesc['lat1'], griddesc['lon1'] = PS.rotatecoords(griddesc['known_lat'],griddesc['known_lon'], \
griddesc['lat0'], griddesc['lon0'], griddesc['stdlon'], -1)
else :
ERROR("Unknown grid type")
sys.exit()
#
# Verify the region chosen
#
nbt, nbj_g, nbi_g = geo.variables["XLONG_M"].shape
ndim = len(geo.variables["XLONG_M"].shape)
if ndim == 3 :
nbt, nbj_g, nbi_g = geo.variables["XLONG_M"].shape
elif ndim == 2 :
nbj_g, nbi_g = geo.variables["XLONG_M"].shape
else :
ERROR("Incompatible shape of coordinate variables")
sys.exit()
if ni < 0 and nj < 0 :
istart = 0
ni = nbi_g
......@@ -136,15 +161,19 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
#
# Extract grid
#
lon_full=np.copy(geo.variables["XLONG_M"][0,jstart:jstart+nj,istart:istart+ni])
lat_full=np.copy(geo.variables["XLAT_M"][0,jstart:jstart+nj,istart:istart+ni])
if ndim == 3 :
lon_full=np.copy(geo.variables["XLONG_M"][0,jstart:jstart+nj,istart:istart+ni])
lat_full=np.copy(geo.variables["XLAT_M"][0,jstart:jstart+nj,istart:istart+ni])
elif ndim == 2 :
lon_full=np.copy(geo.variables["XLONG_M"][jstart:jstart+nj,istart:istart+ni])
lat_full=np.copy(geo.variables["XLAT_M"][jstart:jstart+nj,istart:istart+ni])
#
# Compute resolutions
#
res_lon = np.mean(np.gradient(lon_full, axis=1))
res_lat = np.mean(np.gradient(lat_full, axis=0))
#
elif griddesc['type'] == "RegLonLat" :
else :
# We have a regulat lat/lon grid
nbi_g = geo.variables["lon"].shape[0]
nbj_g = geo.variables["lat"].shape[0]
......@@ -171,11 +200,8 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
griddesc['inilon'] = np.copy(geo.variables["lon"][0])
lat_full=np.transpose(np.tile(np.copy(geo.variables["lat"][jstart:jstart+nj]),(ni,1)))
griddesc['inilat'] = np.copy(geo.variables["lat"][0])
#
else :
ERROR("Unknown grid type")
sys.exit()
#
#
#
return griddesc, lon_full, lat_full, res_lon, res_lat
#
#
......@@ -183,7 +209,14 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
def getland (geo, ist, ni, jst, nj) :
vn=list(v.name for v in geo.variables.values())
if "LANDMASK" in vn :
land=geo.variables["LANDMASK"][0,jst:jst+nj,ist:ist+ni]
ndim = len(geo.variables["LANDMASK"].shape)
if ndim == 3 :
land=geo.variables["LANDMASK"][0,jst:jst+nj,ist:ist+ni]
elif ndim == 2 :
land=geo.variables["LANDMASK"][jst:jst+nj,ist:ist+ni]
else :
ERROR("LANDMASK does not have correct shape")
sys.exit()
elif "Contfrac" in vn :
land=geo.variables["Contfrac"][jst:jst+nj,ist:ist+ni]
elif "elevation" in vn :
......@@ -254,8 +287,12 @@ class ModelGrid :
#
# Define projection
#
if griddesc['type'] == "RegXY" :
proj=PS.LambertC(griddesc['dx'], griddesc['known_lon'], griddesc['known_lat'], griddesc['truelat1'], griddesc['truelat2'], griddesc['stdlon'])
if griddesc['type'] == "lambert_conformal_conic" :
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'], \
griddesc['dlon'], griddesc['dlat'])
elif griddesc['type'] == "RegLonLat" :
proj=PS.RegLonLat(self.res_lon, self.res_lat, griddesc['inilon'], griddesc['inilat'])
else :
......@@ -266,7 +303,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.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)
#
geo.close()
#
......
#
# Produce a lambert Conformal grid as available in WRF
#
import numpy as np
from numba import jit
#
import __main__
import configparser
config=configparser.ConfigParser()
config.read("run.def")
config.read(__main__.ConfigFile)
#
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
rad_per_deg = np.pi/180.
deg_per_rad = 180./np.pi
#
##################################################################
#
# Lambert Conformal projection
#
##################################################################
@jit(nopython = True)
def cone(truelat1, truelat2) :
if (np.abs(truelat1-truelat2) > 0.1) :
......@@ -58,7 +61,6 @@ def sub_LC_ijll( i, j,chi1, chi2, hemi, polei, polej, rebydx, stdlon, cone):
if (lon > +180.) : lon = lon - 360.
if (lon < -180.) : lon = lon + 360.
return np.array([lon, lat])
def LC_llij(polyll, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone) :
polyij=[list(sub_LC_llij(lon, lat, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone)) for lon, lat in polyll]
......@@ -85,8 +87,6 @@ def sub_LC_llij(lon, lat, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon
j = int(np.rint(hemi * j))
return np.array([i,j])
class LambertC:
def __init__(self, dx, known_lon, known_lat, truelat1, truelat2, stdlon):
self.dx = dx
......@@ -128,10 +128,74 @@ class LambertC:
polyij = LC_llij(np.array(polyll), self.hemi, self.truelat1, self.truelat2, self.polei, self.polej, self.rebydx, self.stdlon, self.cone)
return polyij
##################################################################
#
# Cassini or rotated pole projection
#
##################################################################
def rotatecoords(ilat, ilon, lat_np, lon_np, lon_0, direction) :
#
#
phi_np = np.deg2rad(lat_np)
lam_np = np.deg2rad(lon_np)
lam_0 = np.deg2rad(lon_0)
rlat = np.deg2rad(ilat)
rlon = np.deg2rad(ilon)
#
if direction < 0 :
dlam = np.pi - lam_0
else :
dlam = lam_np
#
sinphi = np.cos(phi_np)*np.cos(rlat)*np.cos(rlon-dlam) + np.sin(phi_np)*np.sin(rlat)
cosphi = np.sqrt(1.-sinphi*sinphi)
coslam = np.sin(phi_np)*np.cos(rlat)*np.cos(rlon-dlam) - np.cos(phi_np)*np.sin(rlat)
sinlam = np.cos(rlat)*np.sin(rlon-dlam)
if cosphi != 0. :
coslam = coslam/cosphi
sinlam = sinlam/cosphi
olat = np.rad2deg(np.arcsin(sinphi))
olon = np.rad2deg(np.arctan2(sinlam,coslam)-dlam-lam_0+lam_np)
#
if olon < -180. :
olon = olon + 360.
if olon > 180. :
olon = olon - 360.
#
return olat, olon
#
class Cassini :
def __init__ (self, lon0, lat0, stdlon, lon1, lat1, dlon, dlat) :
self.lonrot = lon0
self.latrot = lat0
self.stdlon = stdlon
# Set up the regular lon/lat projection needed.
self.cyl = RegLonLat(dlon, dlat, lon1, lat1)
def ijll(self, polyij) :
compll = self.cyl.ijll(polyij)
if np.abs(self.latrot) != 90. :
polyll = []
for cll in compll:
polyll.append(list(rotatecoords(cll[1], cll[0], self.latrot, self.lonrot, self.stdlon, 1))[::-1])
else :
polyll = compll
return polyll
def llij(self, polyll) :
if np.abs(self.latrot) != 90. :
compll = []
for ll in polyll :
compll.append(list(rotatecoords(ll[1], ll[0], self.latrot, self.lonrot, self.stdlon, -1))[::-1])
else :
compll = polyll
polyij = self.cyl.llij(compll)
return polyij
##################################################################
#
# Regular Lon/Lat projjection
# Regular Lon/Lat projection
#
##################################################################
......@@ -153,7 +217,8 @@ def RLL_llij(polyll, dlon, dlat, ilon, ilat) :
@jit(nopython = True)
def sub_RLL_llij(lon, lat, dlon, dlat, ilon, ilat):
i = int(np.rint((lon-ilon)/dlon)+1)
dlo = lon-ilon
i = int(np.rint((lon-ilon)%360/dlon)+1)
j = int(np.rint((lat-ilat)/dlat)+1)
return np.array([i,j])
......
......@@ -14,9 +14,10 @@ import time
#
# Gert the information from the configuration file.
#
ConfigFile="run.def"
import configparser
config=configparser.ConfigParser()
config.read("run.def")
config.read(ConfigFile)
nbasmax=config.getint("OverAll", "nbasmax")
numop=config.getint("OverAll", "numop", fallback=100)
OutGraphFile=config.get("OverAll","GraphFile")
......
......@@ -12,11 +12,12 @@ from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import time
#
# Gert the information from the configuration file.
# Get the information from the configuration file.
#
ConfigFile="run.def"
import configparser
config=configparser.ConfigParser()
config.read("run.def")
config.read(ConfigFile)
nbasmax=config.getint("OverAll", "nbasmax")
numop=config.getint("OverAll", "numop", fallback=100)
OutGraphFile=config.get("OverAll","GraphFile")
......
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