Multiple Updates:

- Added MERRA data bias corrections
- Added individual Biomass quantity maps
- CSP works again
- Changed the masking for RoofTop (still needs extra work)
- Shapefile location is updated

Author: thushara2020

code/config.py

@@ -67,7 +67,7 @@ def general_settings():
     param["author"] = "Thushara Addanki"  # the name of the person running the script
     param["comment"] = "Potential Analysis"
     param["path_database_windows"] = "..\..\pyGRETA\Database_KS"    # specify the relative (from the runme.py file) or the absolute path to the database folder
-    param["path_database_linux"] = os.path.expanduser("~") + fs + "database_ks"  # specify path to database on linux
+    param["path_database_linux"] = os.path.expanduser("~") + fs + "pygreta" + fs + "0_database"  # specify path to database on linux
 
     if sys.platform.startswith("win"):
         root = param["path_database_windows"]       # use windows location of database folder
@@ -123,10 +123,11 @@ def scope_paths_and_parameters(paths, param, config_file):
                            skip_blank_lines=True, )  # Import parameters from config_files in folder 'configs'
     input_dict = input_df[1].to_dict()  # Convert dataframe to dict with values from the first column
 
-    paths["subregions"] = PathTemp + input_dict["regions"].replace(" ", "")
     param["region_name"] = input_dict["region_name"].replace(" ", "")  # Name tag of the spatial scope
     param["subregions_name"] = input_dict["subregions_name"] .replace(" ", "") # Name tag of the subregions
     param["country_code"] = input_dict["country_code"].replace(" ", "")
+    paths["subregions"] = PathTemp + "gadm40_" + param["country_code"] + "_shp" + fs + input_dict["regions"].replace(
+        " ", "")
     param["year"] = int(input_dict["year"].replace(" ", ""))  # Convert string 'xxxx' to int
     param["technology"] = input_dict["technology"].replace(" ", "").split(',')  # Creat array by comma separated string
     param["gid"] = input_dict["gid"].replace(" ", "")  # Define spatial level according to GID
@@ -148,7 +149,7 @@ def computation_parameters(param):
     :return param: The updated dictionary param.
     :rtype: dict
     """
-    param["nproc"] = 6
+    param["nproc"] = 10
     param["CPU_limit"] = True
     return param
 
@@ -426,27 +427,27 @@ def osm_areas(param):
 def buffers(param):
 
     buffer = {
-        "snow": 4,
-        "water": 1,
-        "wetland": 1,
+        "snow": 4, #Landuse
+        "water": 1, #Landuse
+        "wetland": 1, #Landuse
 
-        "protected_areas_pv": 1,
-        "protected_areas_windon": 2,
-        "protected_areas_windoff": 4,
+        "protected_areas_pv": 1, #ProtectedAreas
+        "protected_areas_windon": 2, #ProtectedAreas
+        "protected_areas_windoff": 4, #ProtectedAreas
 
-        "airport_windon": 16,
-        "boarder": 2,
+        "airport_windon": 16, #Airports
+        "boarder": 2, #gadm
 
-        "commercial_windon" : 2,
-        "industrial_windon" : 1,
-        "mining" : 1,
-        "military_windon" : 2,
-        "park_pv" : 1,
-        "park_windon" : 3,
-        "recreation_windon" : 1,
+        "commercial_windon" : 2, #osm
+        "industrial_windon" : 1, #osm
+        "mining" : 1, #osm
+        "military_windon" : 2, #osm
+        "park_pv" : 1, #osm
+        "park_windon" : 3, #osm
+        "recreation_windon" : 1, #osm
 
-        "settlement_pv": 1,
-        "settlement_windon" : 4,
+        "settlement_pv": 1, #WSF
+        "settlement_windon" : 4, #WSF
 
         "hydrolakes" : 1,
         "hydrorivers_pv" : 1,
@@ -649,7 +650,7 @@ def csp_parameters(param):
     }
     csp["mask"] = {
         "slope": 20,
-        "lu_suitability": np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0]),
+        "lu_suitability": np.array([1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,1,1,1,1,1,0,0,0,0,1,1,1,0,0]),
         "pa_suitability": np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
     }
     csp["weight"] = {
@@ -992,7 +993,7 @@ def output_folders(paths, param):
 
     # Main output folder
     # paths["region"] = root + "03 Intermediate files" + fs + "Files " + region + fs        # old structure of database
-    paths["region"] = os.path.join(root, "..", "1_results", "Files" + region, "")
+    paths["region"] = os.path.join(root, "..", "1_results", "Files " + region, "")
 
     # Output folder for weather data
     paths["weather_data"] = paths["region"] + "Weather data" + fs
@@ -1049,6 +1050,9 @@ def weather_output_paths(paths, param):
     year = str(param["year"])
     paths["W50M"] = paths["weather_data"] + "w50m_" + year + ".mat"
     paths["CLEARNESS"] = paths["weather_data"] + "clearness_" + year + ".mat"
+    paths["SWGDN"] = paths["weather_data"] + "SWGDN_" + year + ".mat"
+    paths["SWTDN"] = paths["weather_data"] + "SWTDN_" + year + ".mat"
+    paths["SWGDN_real"] = paths["weather_data"] + "SWGDN_real_" + year + ".mat"
     paths["T2M"] = paths["weather_data"] + "t2m_" + year + ".mat"
     paths["W50M_offshore"] = paths["weather_data"] + "w50m_offshore" + year + ".mat"
 
@@ -1101,6 +1105,7 @@ def local_maps_paths(paths, param):
     paths["AREA"] = PathTemp + "_Area.mat"  # Area per pixel in m²
     paths["AREA_offshore"] = PathTemp + "_Area_offshore.mat"  # Area per pixel in m²
     paths["TOPO"] = PathTemp + "_Topography.mat"  # Topography
+    paths["TOPO_low"] = PathTemp + "_Topography_low.mat"
     paths["SLOPE"] = PathTemp + "_Slope.mat"  # Slope
     paths["BATH"] = PathTemp + "_Bathymetry.mat"  # Bathymetry
 
@@ -1283,6 +1288,9 @@ def potential_output_paths(paths, param, tech):
 
     # File name for Biomass
     if tech in ["Biomass"]:
+        paths[tech]["BIOMASS_CROPS"] = PathTemp + "_Biomass_Crops.mat"
+        paths[tech]["BIOMASS_WOOD"] = PathTemp + "_Biomass_Wood.mat"
+        paths[tech]["BIOMASS_MANURE"] = PathTemp + "_Biomass_Manure.mat"
         paths[tech]["BIOMASS_ENERGY"] = PathTemp + "_Biomass_Energy.mat"
         paths[tech]["BIOMASS_CO2"] = PathTemp + "_Biomass_CO2.mat"
     # File name for all other tech
@@ -1339,7 +1347,7 @@ def regional_analysis_output_paths(paths, param, tech):
     paths[tech]["Region_Stats"] = PathTemp + "_Region_stats_" + year + ".csv"
 
     if tech != "Biomass":
-        paths[tech]["Locations"] = PathTemp + "_Locations.shp"
+        paths[tech]["Locations"] = PathTemp + "_Locations_" + year + ".shp"
         paths[tech]["TS"] = PathTemp + "_TS_" + year + ".csv"
         paths[tech]["Sorted_FLH"] = PathTemp + "_sorted_FLH_sampled_" + year + ".mat"
         paths[tech]["Regression_coefficients"] = paths["regression_out"] + subregions + "_" + tech + "_reg_coefficients_"

code/initialization.py

@@ -3,6 +3,9 @@
 from lib.log import logger
 import geopandas as gpd
 import numpy as np
+import urllib.request
+import os
+import zipfile
 import shapely
 # import matplotlib.pyplot as plt
 
@@ -81,7 +84,7 @@ def read_EEZ(paths, param, scope_shp):
     eez_shp = gpd.read_file(paths["EEZ_global"])    # Import all exclusive economic zones worldwide
     eez_shp = eez_shp.to_crs(epsg=4326)             # In case it is not already in this format
 
-    countries = scope_shp['GID_0'].drop_duplicates()        # Obtain the countries out of the regions shapefile
+    countries = scope_shp['ID_0'].drop_duplicates()        # Obtain the countries out of the regions shapefile
     eez_shp = eez_shp[eez_shp['ISO_Ter1'].isin(countries)].reset_index()    # Remain only eez areas of countries involved
 
     param["offshore_scope"] = sf.define_spatial_scope(eez_shp)

code/lib/input_maps.py

@@ -52,15 +52,15 @@ def generate_maps_for_scope(paths, param, multiprocessing):
     """
 
     generate_land(paths, param)  # Subregions
-    generate_sea(paths, param)  # Sea
+    # generate_sea(paths, param)  # Sea
+    generate_topography(paths, param)
     generate_weather_files(paths, param)  # MERRA Weather data for land
-    generate_weather_offshore_files(paths, param)  # MERRA Weather data for offshore # comment out if no offshore
+    # generate_weather_offshore_files(paths, param)  # MERRA Weather data for offshore # comment out if no offshore
 
     if multiprocessing:
         processes = []
 
         processes.append(mp.Process(target=generate_area, args=(paths, param)))
-        processes.append(mp.Process(target=generate_topography, args=(paths, param)))
         processes.append(mp.Process(target=generate_landuse, args=(paths, param)))
         processes.append(mp.Process(target=generate_protected_areas, args=(paths, param)))
         processes.append(mp.Process(target=generate_airports, args=(paths, param)))
@@ -244,8 +244,30 @@ def generate_weather_files(paths, param): #ToDo: Multiprocessing?
 
         # Create the overall wind speed
         W50M = abs(U50M + (1j * V50M))
+        hdf5storage.writes({"SWTDN": SWTDN}, paths["SWTDN"], store_python_metadata=True,
+                           matlab_compatible=True)
+        hdf5storage.writes({"SWGDN": SWGDN}, paths["SWGDN"], store_python_metadata=True,
+                           matlab_compatible=True)
         # Calculate the clearness index
-        CLEARNESS = np.divide(SWGDN, SWTDN, out=np.zeros_like(SWGDN), where=SWTDN != 0)
+        # MERRA bias corrections
+        TOPO_low = hdf5storage.read("TOPO_low", paths["TOPO_low"]).astype(float)
+        TOPO_low = TOPO_low / 1000 #in km
+        SWGDN = hdf5storage.read("SWGDN", paths["SWGDN"]).astype(float)
+        SWTDN = hdf5storage.read("SWTDN", paths["SWTDN"]).astype(float)
+        numfactor1 = 1 + 0.087 * np.multiply(TOPO_low, TOPO_low) - 0.065 * TOPO_low - 0.51
+        # numfactor2 = np.multiply(SWGDN, SWTDN)
+        num = np.zeros([param["m_low"],param["n_low"],8760])
+        for hour in range(8760):
+            num[:,:,hour] = np.multiply(numfactor1, SWGDN[:,:,hour])
+        den = -0.82 * SWGDN + SWTDN
+        CLEARNESS = np.divide(num, den, out=np.zeros_like(num), where=den != 0)
+        for i in range(param["m_low"]):
+            for j in range(param["n_low"]):
+                for k in range(8760):
+                    if CLEARNESS[i,j,k] > 0.7:
+                        CLEARNESS[i,j,k] = SWGDN[i,j,k] / SWTDN[i,j,k]
+        # Clearness without corrections
+        # CLEARNESS = np.divide(SWGDN, SWTDN, out=np.zeros_like(SWGDN), where=SWTDN != 0)
 
         sys.stdout.write("\n")
         logger.info("Writing Files: T2M, W50M, CLEARNESS")
@@ -640,55 +662,27 @@ def generate_topography(paths, param):
         logger.info('Skip')    # Skip generation if files are already there
 
     else:
+    # if 1==1:
         logger.info("Start")
         Crd_all = param["Crd_all"]
-        # Ind = sf.ind_global(Crd_all, param["res_topography"])[0]
         Ind = sf.ind_global(Crd_all, param["res_desired"])[0]
         GeoRef = param["GeoRef"]
-        # Topo = np.zeros((int(180 / param["res_topography"][0]), int(360 / param["res_topography"][1])))
-        # tile_extents = np.zeros((24, 4), dtype=int)
-        # i = 1
-        # j = 1
-        # for letter in ul.char_range("A", "X"):
-        #     north = (i - 1) * 45 / param["res_topography"][0] + 1
-        #     east = j * 60 / param["res_topography"][1]
-        #     south = i * 45 / param["res_topography"][0]
-        #     west = (j - 1) * 60 / param["res_topography"][1] + 1
-        #     tile_extents[ord(letter) - ord("A"), :] = [north, east, south, west]
-        #     j = j + 1
-        #     if j == 7:
-        #         i = i + 1
-        #         j = 1
-        # n_min = (Ind[0] // (45 * 240)) * 45 / param["res_topography"][0] + 1
-        # e_max = (Ind[1] // (60 * 240) + 1) * 60 / param["res_topography"][1]
-        # s_max = (Ind[2] // (45 * 240) + 1) * 45 / param["res_topography"][0]
-        # w_min = (Ind[3] // (60 * 240)) * 60 / param["res_topography"][1] + 1
-        #
-        # need = np.logical_and(
-        #     (np.logical_and((tile_extents[:, 0] >= n_min), (tile_extents[:, 1] <= e_max))),
-        #     np.logical_and((tile_extents[:, 2] <= s_max), (tile_extents[:, 3] >= w_min)),
-        # )
-        #
-        # for letter in ul.char_range("A", "X"):
-        #     index = ord(letter) - ord("A")
-        #     if need[index]:
-        #         with rasterio.open(paths["Topo_tiles"] + "15-" + letter + ".tif") as src:
-        #             tile = src.read()
-        #         Topo[tile_extents[index, 0] - 1 : tile_extents[index, 2], tile_extents[index, 3] - 1 : tile_extents[index, 1]] = tile[0, 0:-1, 0:-1]
-        #
-        # A_TOPO = np.flipud(Topo[Ind[0] - 1 : Ind[2], Ind[3] - 1 : Ind[1]])
-        # A_TOPO = sf.adjust_resolution(A_TOPO, param["res_topography"], param["res_desired"], "mean")
-        # A_TOPO = sf.recalc_topo_resolution(A_TOPO, param["res_topography"], param["res_desired"])
+        m_low = param ["m_low"]
+        n_low = param ["n_low"]
 
         with rasterio.open(paths["Topo_global"]) as src:
             A_TOPO = src.read(1, window=rasterio.windows.Window.from_slices(slice(Ind[0] - 1, Ind[2]),
                                                                             slice(Ind[3] - 1, Ind[1])))
             A_TOPO = np.flipud(A_TOPO)
 
+        A_TOPO_low = np.zeros([m_low,n_low])
+        for i in range(m_low):
+            for j in range(n_low):
+                A_TOPO_low[i,j] = np.sum(A_TOPO[i:i+200,j:j+250])/(200*250)
 
         hdf5storage.writes({"TOPO": A_TOPO}, paths["TOPO"], store_python_metadata=True, matlab_compatible=True)
+        hdf5storage.writes({"TOPO_low": A_TOPO_low}, paths["TOPO_low"], store_python_metadata=True, matlab_compatible=True)
         logger.info("files saved: " + paths["TOPO"])
-        # ul.create_json(paths["TOPO"], param, ["region_name", "Crd_all", "res_topography", "res_desired", "GeoRef"], paths, ["Topo_tiles", "TOPO"])
         ul.create_json(paths["TOPO"], param, ["region_name", "Crd_all", "res_desired", "GeoRef"],
                        paths, ["Topo_global", "TOPO"])
         if param["savetiff_inputmaps"]:
@@ -1408,7 +1402,6 @@ def generate_livestock(paths, param):
     :rtype: None
     """
     logger.info("Start")
-    res_desired = param["res_desired"]
     Crd_all = param["Crd_all"]
     Ind = sf.ind_global(Crd_all, param["res_livestock"])[0]
     GeoRef = param["GeoRef"]
@@ -1421,7 +1414,6 @@ def generate_livestock(paths, param):
                                                                           slice(Ind[3] - 1, Ind[1])))
         A_LS = np.flipud(A_LS)
         A_LS = sf.recalc_livestock_resolution(A_LS, param["res_livestock"], param["res_desired"])
-        # print (np.size(A_LS))
         A_LS[A_LS < 0] = float(0)
         A_LS = np.multiply(A_LS, A_area) / (10 ** 6)
 

code/lib/physical_models.py

@@ -203,8 +203,10 @@ def calc_CF_solar(hour, reg_ind, param, merraData, rasterData, tech):
 
     if tech == "CSP":
         # Wind Speed Corrected at 2m
-        w2m_h = ul.resizem(merraData["W50M"][:, :, hour], m_high, n_high)
-        w2m_h = w2m_h[reg_ind] * rasterData["A_WindSpeed_Corr"][reg_ind]
+        # w2m_h = ul.resizem(merraData["W50M"][:, :, hour], m_high, n_high)
+        w2m_h = merraData["W50M"][:, :, hour]
+        # w2m_h = w2m_h[reg_ind] * rasterData["A_WindSpeed_Corr"][reg_ind]
+        w2m_h = w2m_h[reg_ind_h] * rasterData["A_WindSpeed_Corr"][reg_ind_h]
 
         # Wind Speed cutoff filter:
         windfilter = w2m_h >= csp["Wind_cutoff"]