Merge branch 'master' of gitlab.in2p3.fr:ipsl/lmd/intro/routingpp

StationDiag.py

@@ -3,7 +3,7 @@ import numpy as np
 import datetime
 import sys
 import os
-from inspect import currentframe, getframeinfo
+from inspect import currentframe, getframeinfo, getabsfile
 from netCDF4 import Dataset
 from netCDF4 import stringtochar, chartostring, stringtoarr
 #
@@ -19,6 +19,10 @@ stationsconfig = configparser.ConfigParser()
 stationsconfig.read("stations.def")
 GRDCFile=stationsconfig.get("OverAll","GRDCFile",fallback=None)
 StationList=stationsconfig.get("OverAll","StationList",fallback=[]).split(',')
+if stationsconfig.has_option("OverAll","ExcludedList"):
+    ExcludedList=stationsconfig.get("OverAll","ExcludedList",fallback=[]).split(',')
+else:
+    ExcludedList=[]
 GraphFile=stationsconfig.get("OverAll","GraphFile",fallback=None)
 DiagFile=stationsconfig.get("OverAll","DiagFile",fallback="StationDiag.nc")
 StructureFile = stationsconfig.get("OverAll","StructuresFile",fallback=None)
@@ -34,7 +38,13 @@ if not isinstance(GraphFile, str) or not isinstance(StructureFile, str) :
 if len(StationList) < 1 :
     print("ERROR : no stations provided")
     sys.exit()
+print('Station list:')
 print(StationList)
+if len(ExcludedList) < 1 :
+    print("No excluded stations")
+else: 
+    print('Exclusion list:')
+    print(ExcludedList)
 #
 # Some extra modules so we can start to work
 #
@@ -50,23 +60,53 @@ L = LO.PlacedLocations(GraphFile)
 A=AS.AllStations()
 A.fetchnetcdf(GRDCFile, metaonly=True, region=[L.lonrange,L.latrange])
 #
-# Find the IDs of all the stations listed in stations.def
-#
+# Find the IDs of all the stations and exclusions listed in stations.def
+#
+# Excluded stations/regions
+EX=[]
+for teststn in ExcludedList :
+    sublist = A.FindID(teststn)
+    print("With ", teststn," exclude stations: ",str(sublist))
+    EX.extend(sublist)
+    if not sublist :
+        print("Excluded station selection ", teststn," was already not present in the GRDC file over your regions.")
+EX=list(set(EX)) #avoid duplicates
+#Selected stations/regions
 IDs=[]
+initial_calls=[]
 for teststn in StationList :
-    IDs.append(A.FindID(teststn.strip()))
-    if IDs[-1] < 0 :
-        print("Station ", teststn," could not be found in the GRDC file over your regions.")
+    sublist = A.FindID(teststn)
+    print("With ", teststn," select stations: ",str(sublist))
+    for i in sublist: #avoid duplicates and save the initial call(s)
+        if i in IDs:
+            sublist.remove(i)
+            index=IDs.index(i)
+            initial_calls[index]=initial_calls[index]+', '+teststn
+            continue
+        initial_calls.append(teststn) 
+    IDs.extend(sublist)
+    if not sublist :
+        print("Station selection ", teststn," could not be found in the GRDC file over your regions.")
+if len(IDs)!=len(initial_calls): print('WARNING')
+#Subtraction of ALL recurrences of the excluded stations
+for i in EX:
+    if i in IDs:
+        index=IDs.index(i)
+        del initial_calls[index]
+        del IDs[index]
 # Get additional information from the GRDC file.
 Rivers=[]
 Upstream=[]
+Names = []
 for istn in IDs :
     if istn >= 0 :
         Rivers.append(A.RiverfromID(istn))
         Upstream.append(A.UpstreamfromID(istn))
+        Names.append(A.NamefromID(istn))
     else :
         Rivers.append("-")
         Upstream.append(0)
+        Names.append("_")
 #
 # Find the location within the routing graph of the stations to be diagnosed.
 #
@@ -76,7 +116,7 @@ for i in IDs :
         if i in L.ids :
             Lind.append(L.ids.index(i))
         else :
-            print("Station ", StationList[IDs.index(i)], " (ID =",i,") has not been located in the routing graph.")
+            print("Station ", initial_calls[IDs.index(i)], " (ID =",i,") has not been located in the routing graph.")
             Lind.append(-1)
     else :
         Lind.append(-1) 
@@ -103,9 +143,10 @@ ncvar=[]
 for i,loc in enumerate(Lind) :
     if loc >= 0 :
         vname="ST"+str(IDs[i])
-        print("Creating Variable : ", vname, StationList[i])
+        print("Creating Variable : ", vname, initial_calls[i])
         ncvar.append(nc.createVariable(vname, float, ('y','x')))
-        ncvar[-1].Name=StationList[i]
+        ncvar[-1].InitialCall=initial_calls[i]
+        ncvar[-1].Name=Names[i]
         ncvar[-1].River=Rivers[i]
         ncvar[-1].ObsUpstream=Upstream[i]
         ncvar[-1].Fiindex=L.Fiindex[loc]

Stations/AllStations.py

@@ -2,6 +2,7 @@
 import numpy as np
 import datetime
 import sys
+import re
 from netCDF4 import Dataset
 from netCDF4 import stringtochar, chartostring, stringtoarr
 #
@@ -42,7 +43,19 @@ class AllStations :
     #
     # Inquiry function
     #
-    def FindID(self, name) :
+    def is_integer(self,n):
+        try:
+            float(n)
+        except ValueError:
+            return False
+        else:
+            return float(n).is_integer()
+    def findWholeWord(self,w,s):
+        if re.compile(r'\b({0})\b'.format(w), flags=re.IGNORECASE).search(s):
+            return True
+        else:
+            return False
+    def FindIDfromName(self, name) :
         ID = -1
         for i,stn in enumerate(self.stations) :
             if isinstance(stn.Name, str) :
@@ -50,9 +63,66 @@ class AllStations :
             else :
                 stnname=stn.Name.tolist()
             #
-            if stnname.lower().find(name.lower()) >= 0 :
+            if self.findWholeWord(name.strip(),stnname) :
                 ID = stn.Number
         return ID
+    def FindIDfromCountry(self, river, country) :
+        ID = []
+        for i in self.stations:
+            if isinstance(i.Country,str):
+                stncountry=i.Country
+            else:
+                stncountry=i.Country.tolist()
+            if river!='all':
+                if isinstance(i.River,str):
+                    stnriver=i.River
+                else:
+                    stnriver=i.River.tolist()
+            else:
+                stnriver='all'
+            if self.findWholeWord(country.strip(),stncountry) and self.findWholeWord(river.strip(),stnriver):
+                ID.append(i.Number)
+        return ID
+    def FindIDfromRiver(self, river) :
+        ID = []
+        for i in self.stations:
+            if isinstance(i.River,str):
+                stnriver=i.River
+            else:
+                stnriver=i.River.tolist()
+            if self.findWholeWord(river.strip(),stnriver):
+                ID.append(i.Number)
+        return ID
+    def FindIDfromWMOregion(self, reg) :
+        reg=int(reg)
+        ID = []
+        for i in self.stations:
+            if reg<7:
+                stnreg=i.WMOregion
+            else:
+                stnreg=i.WMOregion*100+i.WMOsubregion
+            if stnreg==reg:
+                ID.append(i.Number)
+        return ID
+    def FindID(self, teststn) :
+        IDs = []
+        if self.is_integer(teststn):
+            IDs.append(int(teststn))
+        elif isinstance(teststn,str):
+            if teststn.find('@')!=-1:
+                selection=teststn.split('@')[0]
+                place=teststn.split('@')[1]
+                if self.is_integer(place):
+                    IDs.extend(self.FindIDfromWMOregion(place.strip()))
+                else:
+                    if len(place)==2 and place.isupper():
+                        IDs.extend(self.FindIDfromCountry(selection.strip(),place.strip()))
+                    else:
+                        IDs.extend(self.FindIDfromRiver(place.strip())) 
+            else: 
+                IDs.append(self.FindIDfromName(teststn.strip()))
+        return IDs
+
     def RiverfromID(self, stnid) :
         allids=[i.Number for i in self.stations]
         if isinstance(self.stations[allids.index(stnid)].River, str) :
@@ -63,6 +133,13 @@ class AllStations :
     def UpstreamfromID(self, stnid) :
         allids=[i.Number for i in self.stations]
         return self.stations[allids.index(stnid)].Area
+    def NamefromID(self, stnid) :
+        allids=[i.Number for i in self.stations]
+        if isinstance(self.stations[allids.index(stnid)].Name, str) :
+            name = self.stations[allids.index(stnid)].Name
+        else :
+            name = self.stations[allids.index(stnid)].Name.tolist()
+        return name
     #
     # Function to dump all stations into a NetCDF file
     #
@@ -479,11 +556,12 @@ def stnnbchecker(C) :
         duploc = np.argmin(np.abs(stnnb-dup))
         dreg = int(stnnb[duploc]/100000)
         dsubreg = int((stnnb[duploc]-dreg*100000)/1000)
-        regstnnb = stnnb[(stnnb >= dup) & (stnnb < dreg*100000+(dsubreg+1)*1000)]
+        regstnnb = stnnb[(stnnb >= dup-1000) & (stnnb < dreg*100000+(dsubreg+1)*1000)]
         if regstnnb.shape[0] < 1 :
             t=dreg*100000+(subreg[i]+1)*1000-1
         else :
-            t= np.max(regstnnb)+1
+            regnum = np.min(regstnnb)
+            t= regnum - 1
         if not any(stnnb == t) :
             stnnb[duploc] = t
             C.stations[duploc].Number = t
@@ -498,12 +576,14 @@ def stnnbchecker(C) :
     if np.count_nonzero(stnnb < 0) > 999 :
         subreg[:] = subreg[:]-1
     mask = np.where(stnnb < 0)
-    for i in np.nditer(mask) :
-        regstnnb=stnnb[(stnnb >= reg[i]*100000+subreg[i]*1000) & (stnnb < (reg[i]+1)*100000)]
+    if len(mask[0])>0:
+      for i in np.nditer(mask) :
+        regstnnb = stnnb[(stnnb >= dup-1000) & (stnnb < dreg*100000+(dsubreg+1)*1000)]
         if regstnnb.shape[0] < 1 :
             t = reg[i]*100000+subreg[i]*1000
         else :
-            t = np.max(regstnnb)+1
+            regnum = np.min(regstnnb)
+            t= regnum - 1
         if not any(stnnb == t) :
             stnnb[i] = t
             C.stations[i].Number = stnnb[i]

Stations/MergeStations.py

@@ -46,15 +46,15 @@ C=AS.stnnbchecker(C)
 #
 print("Number of stations in A = ", A.nbstations," in B = ", B.nbstations, " Merged set has : ", C.nbstations)
 print("Time range of input : A=", A.FirstDate(), "-", A.LastDate()," B=", B.FirstDate(), "-", B.LastDate())
-nbmonth = int((C.LastDate()-C.FirstDate()).days/365)*12+12
+nbmonth = int((C.LastDate().year-C.FirstDate().year)+1)*12
 print("Data coverd by merged dataset :", C.FirstDate(), " to ", C.LastDate(), ' or ', nbmonth, ' month.')
 AllDates=[]
 y=C.FirstDate().year
 m=1
 d=15
 for i in range(nbmonth):
-    dm=i%12
-    dy=int((i-dm)/12)+1
+    dm = i%12
+    dy = i//12
     AllDates.append(datetime.datetime(y+dy, m+dm, d))
 #
 C.dumpnetcdf(AllDates, OutFile, "Merge of "+FileA+" and "+FileB)

Stations/WMOref.py

@@ -36,8 +36,8 @@ class WMOREG:
         """
         # Point is (x,y)-(lon,lat) 
         p1 = Point(lon, lat)
-        out = np.array([p1.within(poly) for poly in self.polygonWMO])
-        try : 
+        out = np.array([poly.contains(p1) for poly in self.polygonWMO])
+        try :
             index = np.where(out)[0][0]
             stReg, stSubreg = self.attributesPolygonWMO[index]
             if output=="num":
@@ -46,5 +46,19 @@ class WMOREG:
                 return [stReg, stSubreg]
         except:
             print("lon: {0}; lat: {1}".format(lon,lat))
-            print("No valid WMO region has been found, verify the lon lat")
+            print("Inexact location : finding the closest subregion.")
+            out = np.array([self.get_distance(poly,p1) for poly in self.polygonWMO])
+            index = np.argmin(out)
+            stReg, stSubreg = self.attributesPolygonWMO[index]
+            if output=="num":
+                return stReg*100+stSubreg
+            elif output == "list":
+                return [stReg, stSubreg]
+            
+    def get_distance(self, poly,p1):
+        try:
+           d = poly.exterior.distance(p1)
+        except:
+           d = np.min([p.exterior.distance(p1) for p in poly])
+        return d