update regionalised_dailyT file version

.github/workflows/build.yml

@@ -81,7 +81,7 @@ jobs:
           cat .gitignore_extra >> .gitignore
           pip install numpy==${{matrix.numpy-version}} build==1.2.2.post1 ninja meson
           cd src/pyraingen/fortran_daily
-          python -m numpy.f2py -c regionalised_dailyT.for -m regionalised_dailyT -I $PWD --fcompiler=gfortran --compiler=msvc --backend distutils
+          python -m numpy.f2py -c regionalised_dailyT4.for -m regionalised_dailyT4 -I $PWD --fcompiler=gfortran --compiler=msvc --backend distutils
           Get-ChildItem -Force -R
           cd ../../..
           python -m build
@@ -92,7 +92,7 @@ jobs:
           cat .gitignore_extra >> .gitignore
           pip install numpy==${{matrix.numpy-version}} ninja meson setuptools==59.8.0 build # wheel
           cd src/pyraingen/fortran_daily
-          python -m numpy.f2py -c regionalised_dailyT.for -m regionalised_dailyT --fcompiler=gnu95 --compiler=mingw32
+          python -m numpy.f2py -c regionalised_dailyT4.for -m regionalised_dailyT4 --fcompiler=gnu95 --compiler=mingw32
           Get-ChildItem -Force -R
           cd ../../..
           python -m build

docs/example.ipynb

@@ -63,7 +63,7 @@
     "\n",
     "Use the inbuilt function `pyraingen.get_for_path()` to get the path to the fortran script and change, or point, to that directory to compile. Use the following commands in the cmd prompt.\n",
     "\n",
-    "`python -m numpy.f2py -c regionalised_daily.for -m regionalised_daily --fcompiler=gnu95 --compiler=mingw32`\n",
+    "`python -m numpy.f2py -c regionalised_dailyT4.for -m regionalised_dailyT4 --fcompiler=gnu95 --compiler=mingw32`\n",
     "\n",
     "\n",
     "Point to the location(s) where you have stored the example data. Please copy the `stn_record.dat` from the example daily data to the working directory where the `regionaliseddailysim` function will be executed from."

src/pyraingen/fortran_daily/CMakeLists.txt

@@ -6,8 +6,8 @@ find_package(
   REQUIRED
 )
 
-set(f2py_module_name "regionalised_dailyT")
-set(fortran_src_file "${CMAKE_CURRENT_SOURCE_DIR}/regionalised_dailyT.for")
+set(f2py_module_name "regionalised_dailyT4")
+set(fortran_src_file "${CMAKE_CURRENT_SOURCE_DIR}/regionalised_dailyT4.for")
 
 
 if (WIN32)
@@ -19,7 +19,7 @@ if (WIN32)
   #  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
   #  )
   FILE(GLOB generated_module_file  *.pyd)
-  FILE(GLOB generated_dll  regionalised_dailyT/.libs/*.dll)
+  FILE(GLOB generated_dll  regionalised_dailyT4/.libs/*.dll)
   FILE(GLOB out_files Release/*.o)
   FILE(GLOB module_files  src.win-amd64-*/*)
   # find_package(F2PY REQUIRED)

src/pyraingen/get_for_path.py

@@ -8,6 +8,6 @@ def get_for_path():
     pathlib.PosixPath
         Path to file.
     """
-    with resources.path("pyraingen.fortran_daily", "regionalised_daily.for") as f:
+    with resources.path("pyraingen.fortran_daily", "regionalised_dailyT4.for") as f:
         data_file_path = f
     return str(data_file_path)

src/pyraingen/regionaliseddailysim.py

@@ -18,7 +18,8 @@ def regionaliseddailysim(nyears, startyear, nsim,
                             output_path_nc='daily.nc',
                             output_stats='stat_.out', output_val='rev_.out',
                             cutoff=0.30, wind=15, nstation=5, 
-                            nlon=3, lag=1, iamt=1, ival=0, irf=1, interactive=True):
+                            nlon=3, lag=1, iamt=1, ival=0, irf=1, rm=1.0,
+                            getstations=True):
     """Front end to the regionalised daily rainfall generator.
     For a detailed description of the approach refer to: 
     Mehrotra et al (2012) "Continuous Rainfall simulation: 2.
@@ -90,10 +91,13 @@ def regionaliseddailysim(nyears, startyear, nsim,
         Default is 0.
     irf : int
         Leave at default.
-        Default is 1. 
-    interactive : bool
-        If True, user will be asked to drop nearby stations interactively.
-        Default is True.
+        Default is 1.
+    rm : float
+        Rainfall multiplier to scale mean annual rainfall.
+        Default is 1.0
+    getstations : boolean
+        Get nearby stations or use existing "nearby_station_details.out"
+        Default is True        
 
     Returns
     ----------   
@@ -129,7 +133,7 @@ def regionaliseddailysim(nyears, startyear, nsim,
     with open("data_r.dat",'r') as file:
         data_r = file.readlines()
 
-    data_r[2] = ' ' + rain + iband + nstation + nsim + nlon + lag + ng + nsgtart + iamt + ival + irf +'\n'
+    data_r[2] = ' ' + rain + iband + nstations + nsim + nlon + lag + ng + nsgtart + iamt + ival + irf +'\n'
     data_r[10] = ' ' + idx + lat + lon + elev + dcoast + anrf + temp + '\n'
     data_r[12] = data_path +'\n'
     data_r[14] = f"mmm_{targetidx}.out" +'\n'
@@ -139,44 +143,96 @@ def regionaliseddailysim(nyears, startyear, nsim,
     with open("data_r.dat",'w') as file:
         file.writelines(data_r)
 
+
+    if getstations:
+        # Check if file exists
+        if os.path.exists("nearby_station_details.out"):
+            os.remove("nearby_station_details.out")
+
+        # Get Data
+        if pathStnData == None:
+            with resources.path("pyraingen.data", "stn_record.csv") as f:
+                pathStnData = str(f)
+        if pathModelCoeffs == None:
+            with resources.path("pyraingen.data", "daily_logreg_coefs.csv") as f:
+                pathModelCoeffs = str(f)
+
+        # Get Nearby Stations
+        param = {}
+        param['nNearStns']       = nstation
+        param['pathStnData']     = pathStnData
+        param['pathModelCoeffs'] = pathModelCoeffs
+        param['pathDailyData']   = data_path
+
+        target = {}
+        target['index']           = targetidx             
+        target['lat']             = targetlat               
+        target['lon']             = targetlon               
+        target['elevation']       = targetelev              
+        target['distToCoast']     = targetdcoast       
+        target['annualRainDepth'] = targetanrf
+        target['temp']            = targettemp    
+
+        print("\nFinding nearby stations...\n")
+        station(param, target, nAttributes=33,
+            fout=f'nearby_station_details.out'
+        )
+
+    # Check if file exists
+    # if os.path.exists("drop.out"):
+        # os.remove("drop.out")
+
     # Begin simulation
-    idrop =-1
-    kk = 0
-    break_out_flag = False
+    regionalised_dailyT4.regionalised_daily(idrop=0)
+
+    # idrop = 0
+    # kk = 0
+    # break_out_flag = False
 
-    if interactive:
-        while True:
-            regionalised_dailyT.regionalised_daily(idrop=idrop)
-            with open('nearby_station_details.out') as file:
-                print(file.read())
-            nearby = np.genfromtxt('nearby_station_details.out', 
-                                    skip_header=6, invalid_raise=False) 
-            while True:
-                    try:
-                        idrop = int(input("Do you want to drop any nearby station?\nIf yes, enter station sr no otherwise enter zero: " ))
-                        if idrop < 0:
-                            print("Sorry, no numbers below zero.") 
-                        elif idrop > 0:
-                            if idrop > 5:
-                                print("Sorry, please choose one of the listed nearby station sr nos.")
-                            else:
-                                kk += 1
-                                with open ("drop.out", 'a') as file:
-                                    file.write(str(int(kk)).rjust(12)
-                                                + str(int(nearby[idrop-1][1])).rjust(12)+'\n')
-                                idrop = kk
-                                break
-                        else:
-                            regionalised_dailyT.regionalised_daily(idrop=idrop)
-                            break_out_flag = True
-                            break
-                    except ValueError:
-                        print("Sorry, invalid input. Please makesure input is integer only.")
-            if break_out_flag:
-                break
+    # while True:
+        # regionalised_dailyT.regionalised_daily(idrop=idrop)
+        # with open('nearby_station_details.out') as file:
+            # print(file.read())
+        # nearby = np.genfromtxt('nearby_station_details.out', 
+                                # skip_header=6, invalid_raise=False) 
+        # while True:
+                # try:
+                    # idrop = int(input("Do you want to drop any nearby station?\nIf yes, enter station sr no otherwise enter zero: " ))
+                    # if idrop < 0:
+                        # print("Sorry, no numbers below zero.") 
+                    # elif idrop > 0:
+                        # if idrop > 5:
+                            # print("Sorry, please choose one of the listed nearby station sr nos.")
+                        # else:
+                            # kk += 1
+                            # with open ("drop.out", 'a') as file:
+                                # file.write(str(int(kk)).rjust(12)
+                                            # + str(int(nearby[idrop-1][1])).rjust(12)+'\n')
+                            # idrop = kk
+                            # break
+                    # else:
+                        # regionalised_dailyT.regionalised_daily(idrop=idrop)
+                        # break_out_flag = True
+                        # break
+                # except ValueError:
+                    # print("Sorry, invalid input. Please makesure input is integer only.")
+
+        # if break_out_flag:
+            # break
+
+
+    convertdailync(f"mmm_{targetidx}.out", output_path_nc, startyear, nyears, nsim, missingDay=-999.9)
+    #remove textfile
+    if os.path.exists(f"mmm_{targetidx}.out"):
+        os.remove(f"mmm_{targetidx}.out")
 
-    if netcdf == True:
-        convertdailync(output_path_txt, output_path_nc, startyear, nyears, nsim, missingDay=-999.9)
-        #remove textfile
-        if os.path.exists(output_path_txt):
-            os.remove(output_path_txt)
+    # Bias Correct and/or Scale Rainfall
+    daily_rain = xr.open_dataset(output_path_nc)
+    daily_rain['day'] = pd.to_datetime(daily_rain['day'], unit='D', origin='julian')
+    smanrf = daily_rain['rainfall'].resample(day="A").sum().mean()
+    daily_rain.close()
+    del daily_rain
+    daily_rain_bc = xr.open_dataset(output_path_nc)    
+    daily_rain_bc['rainfall'] *= (targetanrf * rm)/smanrf
+    daily_rain_bc.close()
+    daily_rain_bc.to_netcdf(output_path_nc)