Include Models Folder

Collect_Tools.py

@@ -34,6 +34,7 @@ def Newest():
     wa_folder_Product = os.path.join(home_folder,'Products') 
     wa_folder_Generator = os.path.join(home_folder,'Generator') 
     wa_folder_Functions = os.path.join(home_folder,'Functions') 
+    wa_folder_Models = os.path.join(home_folder,'Models') 	
     wa_folder_Sheets = os.path.join(home_folder,'Sheets') 
 
 
@@ -43,6 +44,7 @@ def Newest():
     wa_master_folder_Collect = os.path.join(home_folder, 'wa-master','Collect') 
     wa_master_folder_Product = os.path.join(home_folder, 'wa-master','Products') 
     wa_master_folder_Sheets = os.path.join(home_folder, 'wa-master','Sheets') 
+    wa_master_folder_Models = os.path.join(home_folder, 'wa-master','Models') 	
     wa_master_folder_Home = os.path.join(home_folder, 'wa-master') 
 
     shutil.rmtree(wa_folder_General)
@@ -51,13 +53,15 @@ def Newest():
     shutil.rmtree(wa_folder_Sheets)
     shutil.rmtree(wa_folder_Generator)
     shutil.rmtree(wa_folder_Functions)
+    shutil.rmtree(wa_folder_Models)
 
     shutil.copytree(wa_master_folder_General, wa_folder_General)
     shutil.copytree(wa_master_folder_Collect, wa_folder_Collect)
     shutil.copytree(wa_master_folder_Product, wa_folder_Product)
     shutil.copytree(wa_master_folder_Sheets, wa_folder_Sheets)
     shutil.copytree(wa_master_folder_Generator, wa_folder_Generator)
     shutil.copytree(wa_master_folder_Functions, wa_folder_Functions)
+    shutil.copytree(wa_master_folder_Models, wa_folder_Models)
 
     shutil.rmtree(wa_master_folder_Home)				
     os.remove(file_nametext)

Functions/Five/Irrigation.py

@@ -35,10 +35,10 @@ def Calc_Supply_Budyko(Discharge_dict, Name_NC_Rivers, Name_NC_ET, Name_NC_ETref
     return(DataCube_surface_withdrawal_m3, DataCube_ETblue_m3)
 
 def Add_irrigation(Discharge_dict, River_dict, Name_NC_Rivers, Name_NC_Supply, Name_NC_ETblue, Name_NC_Basin, Startdate, Enddate, Example_dataset):
-
+    #Discharge_dict_CR2, River_dict_CR2, Name_NC_Rivers_CR,  Name_NC_Supply, Name_NC_ETblue, Name_NC_Basin_CR, Startdate, Enddate, Example_dataset
     import copy
     import numpy as np
-
+    import sys
     import wa.General.raster_conversions as RC
 
     # Copy dicts as starting adding reservoir 
@@ -63,20 +63,27 @@ def Add_irrigation(Discharge_dict, River_dict, Name_NC_Rivers, Name_NC_Supply, N
     # find IDs drainage for only the basin
     Basin = RC.Open_nc_array(Name_NC_Basin)
     ID_Rivers_flow = RC.gap_filling(ID_Rivers,NoDataValue = 0.) * Basin
-
-    for i in np.unique(ID_Rivers_flow):
+    Water_Error = 0 
+
+    for i in np.unique(ID_Rivers_flow)[1:]:
         if np.nansum(DataCube_ETblue_m3[:,ID_Rivers_flow == i]) > 0:
+            sys.stdout.write("\r%s Procent of adding irrigation completed with %s x 10^9 m3 Water Error       " %(np.int((i-np.unique(ID_Rivers_flow)[1])/(np.unique(ID_Rivers_flow)[-1]-np.unique(ID_Rivers_flow)[1])*100),Water_Error/1e9))
+            sys.stdout.flush()
+            #print('%s Procent of adding irrigation completed' %np.int(i/np.unique(ID_Rivers_flow)[-1]*100))
             total_surface_withdrawal = np.nansum(DataCube_surface_withdrawal_m3[:,ID_Rivers_flow == i] ,1)
 
             # Find exact area in river directory
             for River_part in River_dict.iteritems():
                 if len(np.argwhere(River_part[1] == i)) > 0:
-                    
+
                     # Find the river part in the dictionery
                     row_discharge = np.argwhere(River_part[1]==i)[0][0]
 
                     # Subtract the withdrawal from that specific riverpart
-                    Discharge_dict_new[River_part[0]][:,0:row_discharge] = Discharge_dict_new[River_part[0]][:,0:row_discharge] - total_surface_withdrawal[:,None] 
+                    Real_Surface_Withdrawal = np.minimum(Discharge_dict_new[River_part[0]][:,row_discharge].flatten(), total_surface_withdrawal[:,None].flatten())
+
+                    Water_Error += np.maximum(np.nansum(total_surface_withdrawal[:,None].flatten() - Real_Surface_Withdrawal),0)
+                    Discharge_dict_new[River_part[0]][:,0:row_discharge] = Discharge_dict_new[River_part[0]][:,0:row_discharge] - Real_Surface_Withdrawal[:,None] 
 
                     # Subtract the withdrawal from the part downstream of the riverpart within the same dictionary
                     Discharge_dict_new[River_part[0]][np.logical_and(Discharge_dict_new[River_part[0]]<=0,Discharge_dict[River_part[0]]>=0)] = 0
@@ -86,10 +93,12 @@ def Add_irrigation(Discharge_dict, River_dict, Name_NC_Rivers, Name_NC_Supply, N
                     # Subtract the withdrawal from all the other downstream dictionaries
                     while len(River_dict) > times:
                         for River_part_downstream in River_dict.iteritems():
-                            if River_dict[River_part[0]][-1] == End_river:  
-                                print River_part_downstream
-                                Discharge_dict_new[River_part_downstream[0]][:,1:] = Discharge_dict_new[River_part_downstream[0]][:,1:] - total_surface_withdrawal[:,None] 
-                                Discharge_dict_new[River_part[0]][np.logical_and(Discharge_dict_new[River_part[0]]<=0,Discharge_dict[River_part[0]]>=0)] = 0
+                            if River_part_downstream[1][-1] == End_river:  
+
+                                Discharge_dict_new[River_part_downstream[0]][:,:] = Discharge_dict_new[River_part_downstream[0]][:,:] - Real_Surface_Withdrawal[:,None] 
+                                #Discharge_dict_new[River_part_downstream[0]][:,1:] = Discharge_dict_new[River_part_downstream[0]][:,1:] - total_surface_withdrawal[:,None] 
+
+                                Discharge_dict_new[River_part_downstream[0]][np.logical_and(Discharge_dict_new[River_part_downstream[0]]<=0,Discharge_dict[River_part_downstream[0]]>=0)] = 0
                                 End_river = River_dict[River_part_downstream[0]][0]
                                 times = 0
                             times += 1

General/data_conversions.py

@@ -251,14 +251,16 @@ def Save_as_NC(namenc, DataCube, Var, Reference_filename,  Startdate = '', Endda
             timeo.standard_name = 'time'
 
         # Create the lon variable
-        lono = nco.createVariable('longitude', 'f4', ('longitude',))
+        lono = nco.createVariable('longitude', 'f8', ('longitude',))
         lono.standard_name = 'longitude'
         lono.units = 'degrees_east'
+        lono.pixel_size = geo_out[1]
 
         # Create the lat variable
-        lato = nco.createVariable('latitude', 'f4', ('latitude',))
+        lato = nco.createVariable('latitude', 'f8', ('latitude',))
         lato.standard_name = 'latitude'
         lato.units = 'degrees_north'
+        lato.pixel_size = geo_out[5]
 
         # Create container variable for CRS: lon/lat WGS84 datum
         crso = nco.createVariable('crs', 'i4')

General/raster_conversions.py

@@ -78,8 +78,8 @@ def Open_nc_info(NC_filename):
     lats = fh.variables['latitude'][:]
     lons = fh.variables['longitude'][:]
 
-    Geo6 = lats[1]-lats[0]
-    Geo2 = lons[1]-lons[0]
+    Geo6 = fh.variables['latitude'].pixel_size
+    Geo2 = fh.variables['longitude'].pixel_size
     Geo4 = np.max(lats) + Geo6/2	
     Geo1 = np.min(lons) - Geo2/2
 
@@ -301,7 +301,7 @@ def reproject_dataset_example(dataset, dataset_example, method=1):
 
     # open dataset that must be transformed 
     try:
-        if dataset.split('.')[-1] == 'tif':
+        if os.path.splitext(dataset)[-1] == '.tif':
             g = gdal.Open(dataset)               					
         else:
             g = dataset    
@@ -313,13 +313,13 @@ def reproject_dataset_example(dataset, dataset_example, method=1):
     if epsg_from == 9001:
         epsg_from = 5070
 
-
     # open dataset that is used for transforming the dataset
     try:
-        if dataset_example.split('.')[-1] == 'tif':
-            gland = gdal.Open(dataset_example)               					
+        if os.path.splitext(dataset_example)[-1] == '.tif':
+            gland = gdal.Open(dataset_example)   	
         else:
             gland = dataset_example      
+
     except:
             gland = dataset_example              
     epsg_to = Get_epsg(gland)	
@@ -495,9 +495,29 @@ def Get3Darray_time_series_monthly(Dir_Basin, Data_Path, Startdate, Enddate, Exa
         file_name_path = os.path.join(Dir_Basin, Data_Path, file_name[0])								
         if Example_data is not None:
             if Date == Dates[0]:
+                if os.path.splitext(Example_data)[-1] == '.tif':
                     geo_out, proj, size_X, size_Y = Open_array_info(Example_data)													
                     dataTot=np.zeros([len(Dates),size_Y,size_X])														
-
+                if os.path.splitext(Example_data)[-1] == '.nc':
+                    geo_out, projection, size_X, size_Y, size_Z, Time = Open_nc_info(Example_data)													
+                    dataTot=np.zeros([len(Dates),size_Y,size_X])	 
+
+                    # Create memory file for reprojection
+                    data = Open_nc_array(Example_data)
+                    driver = gdal.GetDriverByName("MEM")
+                    gland = driver.Create('', size_Y, size_X, 1,
+                                           gdal.GDT_Float32, ['COMPRESS=LZW'])
+                    srse = osr.SpatialReference()
+                    if projection == '' or projection == 4326:
+                        srse.SetWellKnownGeogCS("WGS84")
+                    else:	
+                        srse.SetWellKnownGeogCS(projection)
+                    gland.SetProjection(srse.ExportToWkt())
+                    gland.GetRasterBand(1).SetNoDataValue(-9999)
+                    gland.SetGeoTransform(geo_out)
+                    gland.GetRasterBand(1).WriteArray(data)		
+                    Example_data = gland
+
             dest = reproject_dataset_example(file_name_path, Example_data, method=1)
             Array_one_date = dest.GetRasterBand(1).ReadAsArray() 								
         else: 								

Generator/Sheet5/main.py

@@ -218,19 +218,23 @@ def Calculate(WA_HOME_folder, Basin, P_Product, ET_Product, Inflow_Text_Files, W
 
             # Calculate runoff based on Budyko
             DataCube_Runoff_CR = Five.Budyko.Calc_runoff(Name_NC_ETref_CR, Name_NC_Prec_CR)
-
+
+            # Save the runoff as netcdf
+            DC.Save_as_NC(Name_NC_Runoff_for_Routing_CR, DataCube_Runoff_CR, 'Runoff_CR', Example_dataset, Startdate, Enddate, 'monthly', 0.01)
+            del DataCube_Runoff_CR    
+
     ###################### 5b. Get Runoff from WaterPIX ###########################
     if WaterPIX_filename != "": 
 
-        # Get WaterPIX data
-        WaterPIX_Var = 'TotalRunoff_M'
-        DataCube_Runoff_CR = Five.Read_WaterPIX.Get_Array(WaterPIX_filename, WaterPIX_Var, Example_dataset, Startdate, Enddate)
-                  
-
-    # Save the runoff as netcdf
-    DC.Save_as_NC(Name_NC_Runoff_for_Routing_CR, DataCube_Runoff_CR, 'Runoff_CR', Example_dataset, Startdate, Enddate, 'monthly', 0.01)
-    del DataCube_Runoff_CR    
-
+        if not os.path.exists(Name_NC_Runoff_CR):       
+            # Get WaterPIX data
+            WaterPIX_Var = 'TotalRunoff_M'
+            DataCube_Runoff_CR = Five.Read_WaterPIX.Get_Array(WaterPIX_filename, WaterPIX_Var, Example_dataset, Startdate, Enddate)
+                      
+            # Save the runoff as netcdf
+            DC.Save_as_NC(Name_NC_Runoff_for_Routing_CR, DataCube_Runoff_CR, 'Runoff_CR', Example_dataset, Startdate, Enddate, 'monthly', 0.01)
+            del DataCube_Runoff_CR    
+         
     '''  
     ###################### Calculate Runoff with P min ET ###########################
   
@@ -309,7 +313,7 @@ def Calculate(WA_HOME_folder, Basin, P_Product, ET_Product, Inflow_Text_Files, W
         # Define the 2% highest pixels as rivers
         Rivers = np.zeros([np.size(Routed_Discharge_Ave,0),np.size(Routed_Discharge_Ave,1)])
         Routed_Discharge_Ave[Raster_Basin != 1] = np.nan
-        Routed_Discharge_Ave_number = np.nanpercentile(Routed_Discharge_Ave,90)
+        Routed_Discharge_Ave_number = np.nanpercentile(Routed_Discharge_Ave,99)
         Rivers[Routed_Discharge_Ave > Routed_Discharge_Ave_number] = 1  # if yearly average is larger than 5000km3/month that it is a river
 
         # Save the river file as netcdf file
@@ -410,18 +414,17 @@ def Calculate(WA_HOME_folder, Basin, P_Product, ET_Product, Inflow_Text_Files, W
 
     Name_py_Discharge_dict_CR3 = os.path.join(Dir_Basin,'Simulations','Simulation_%d' %Simulation, 'Sheet_5','Discharge_dict_CR3_simulation%d.npy' %(Simulation))    
     Name_NC_Supply = DC.Create_NC_name('Supply', Simulation, Dir_Basin, 5, info)
-    Name_NC_ETblue =  DC.Create_NC_name('ETblue', Simulation, Dir_Basin, 5, info)
 
     if not os.path.exists(Name_py_Discharge_dict_CR3):
 
         if Supply_method == "Fraction":
+            Name_NC_ETblue =  DC.Create_NC_name('ETblue', Simulation, Dir_Basin, 5, info)
             DataCube_Supply_m3, DataCube_ETblue_m3 = Five.Irrigation.Calc_Supply_Budyko(Discharge_dict_CR2, Name_NC_Rivers_CR, Name_NC_ET_CR, Name_NC_ETref_CR, Name_NC_Prec_CR, Name_NC_Basin_CR, Name_NC_frac_sw_CR, Startdate, Enddate, Example_dataset)
             DC.Save_as_NC(Name_NC_Supply, DataCube_Supply_m3, 'Supply', Example_dataset, Startdate, Enddate, 'monthly')
             DC.Save_as_NC(Name_NC_ETblue, DataCube_ETblue_m3, 'ETblue', Example_dataset, Startdate, Enddate, 'monthly') 
             del DataCube_ETblue_m3, DataCube_Supply_m3
             Discharge_dict_CR3 = Five.Irrigation.Add_irrigation(Discharge_dict_CR2, River_dict_CR2, Name_NC_Rivers_CR,  Name_NC_Supply, Name_NC_ETblue, Name_NC_Basin_CR, Startdate, Enddate, Example_dataset)
-            np.save(Name_py_Discharge_dict_CR3, Discharge_dict_CR3)
-
+            np.save(Name_py_Discharge_dict_CR3, Discharge_dict_CR3)           
 
         if Supply_method == "WaterPIX":
             WaterPIX_Var = 'Supply_M'
@@ -431,6 +434,13 @@ def Calculate(WA_HOME_folder, Basin, P_Product, ET_Product, Inflow_Text_Files, W
             DC.Save_as_NC(Name_NC_Supply, DataCube_Supply_m3, 'Supply', Example_dataset, Startdate, Enddate, 'monthly')
             Discharge_dict_CR3 = Five.Irrigation.Add_irrigation(Discharge_dict_CR2, River_dict_CR2, Name_NC_Rivers_CR,  Name_NC_Supply, Name_NC_Supply, Name_NC_Basin_CR, Startdate, Enddate, Example_dataset)
 
+        if Supply_method == "UserInput":
+
+            # Get the data of Supply defined by user and save as nc
+            DataCube_Supply_m3 = RC.Get3Darray_time_series_monthly(Dir_Basin, WaterPIX_filename, Startdate, Enddate, Example_data = Example_dataset)
+            DC.Save_as_NC(Name_NC_Supply, DataCube_Supply_m3, 'ETref_CR', Example_dataset, Startdate, Enddate, 'monthly', 0.01)
+            Discharge_dict_CR3 = Five.Irrigation.Add_irrigation(Discharge_dict_CR2, River_dict_CR2, Name_NC_Rivers_CR,  Name_NC_Supply, Name_NC_Supply, Name_NC_Basin_CR, Startdate, Enddate, Example_dataset)
+
     else:
         Discharge_dict_CR3 = np.load(Name_py_Discharge_dict_CR3).item() 
 
@@ -456,6 +466,34 @@ def Calculate(WA_HOME_folder, Basin, P_Product, ET_Product, Inflow_Text_Files, W
 
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
     '''
 
     # DEM