Backup of first investigation

Author: stelliom

python/extpar_art_to_buffer.py

@@ -2,6 +2,7 @@
 import logging
 import netCDF4 as nc
 import numpy as np
+from time import perf_counter
 
 from joblib import Parallel, delayed, dump, load
 from tqdm import tqdm
@@ -95,6 +96,81 @@ def get_fraction_per_soil_type(lu):
     return fracs
 
 
+def calculate_soil_fraction_optimized(target_grid, soil_types_raw, nearest_target_cell_to_raw_cells, ncpu=2):
+    """
+    target_grid: target ICON grid
+    soil_types_raw: landuse class for each cell from the HWSD dataset (LU variable)
+    nearest_target_cell_to_raw_cell: indices of the cell from the target ICON grid which is nearest to each cell of the raw grid (from HWSD dataset)
+    """
+    ncells_target = target_grid.lons.size
+    nsoil_types = 13
+
+    soil_ids = np.arange(1, nsoil_types+1)
+    soil_fractions_target = np.zeros((ncells_target, nsoil_types))
+
+    target_cells, n_nearest_raw_cells = np.unique(nearest_target_cell_to_raw_cells, return_counts=True)
+
+    for soil_id in tqdm(soil_ids[:1]):
+
+        target_cells_with_soil_type, n_nearest_raw_cells_with_soil_type = np.unique(np.where(soil_types_raw == soil_id, nearest_target_cell_to_raw_cells, -1), return_counts=True)
+
+        for target_cell_id in np.arange(ncells_target):
+
+            soil_fraction = np.array(n_nearest_raw_cells_with_soil_type[target_cells_with_soil_type == target_cell_id] / n_nearest_raw_cells[target_cells == target_cell_id])
+
+            if len(soil_fraction) != 0:
+                soil_fractions_target[target_cell_id, soil_id - 1] = soil_fraction
+
+    return soil_fractions_target
+
+
+def calculate_soil_fraction_test(tg, lus, idxs, ncpu=2):
+    """
+    lus: LU classes from HWSD data
+    idxs: indices corrsponding to icon grid for each grid in HWSD data
+    tg: ICON grid
+    """
+    soil_types = np.arange(1, 14)
+    fracs = np.zeros((tg.lons.size, soil_types.size))
+    grid_ids, grid_counts = np.unique(idxs, return_counts=True)
+    print("soil_types:", soil_types.shape)
+    print("fracs:", fracs.shape)
+    print("grid_ids:", grid_ids.shape)
+    print("grid_counts:", grid_counts.shape)
+
+    def get_fraction_per_soil_type(lu):
+        print("lus:", lus.shape)
+        print("idxs:", idxs.shape)
+        test=np.where(lus == lu, idxs, -1)
+        print("test:", test.shape)
+        start=perf_counter()
+        grid_class, grid_count = np.unique(np.where(lus == lu, idxs, -1),
+                                           return_counts=True)
+        end=perf_counter()
+        print("grid_class:", grid_class.shape)
+        print("grid_count:", grid_count.shape)
+        print("Sect 1:", end-start)
+        start=perf_counter()
+        for grid_id in grid_class:
+            frac = np.array(grid_count[grid_class == grid_id] /
+                            grid_counts[grid_ids == grid_id])
+            if len(frac) != 0:
+                fracs[grid_id, lu - 1] = frac
+        end=perf_counter()
+        print("Sect 2:", end-start)
+
+    #Parallel(n_jobs=13,
+    #         max_nbytes='100M',
+    #         mmap_mode='w+',
+    #         backend='threading')(delayed(get_fraction_per_soil_type)(lu)
+    #                              for lu in tqdm(soil_types))
+    
+    for lu in tqdm(soil_types[:1]):
+        get_fraction_per_soil_type(lu)
+
+    return fracs
+
+
 # --------------------------------------------------------------------------
 # --------------------------------------------------------------------------
 # initialize logger
@@ -171,6 +247,22 @@ def get_fraction_per_soil_type(lu):
 lons = np.array(tg.lons)
 lats = np.array(tg.lats)
 
+vlons = np.array(tg.vlons)
+vlats = np.array(tg.vlats)
+
+lon_min = np.min(vlons)
+lon_max = np.max(vlons)
+lat_min = np.min(vlats)
+lat_max = np.max(vlats)
+
+lon_mask = (raw_lon >= lon_min) & (raw_lon <= lon_max)
+lat_mask = (raw_lat >= lat_min) & (raw_lat <= lat_max)
+
+raw_lon = raw_lon[lon_mask]
+raw_lat = raw_lat[lat_mask]
+raw_lus = raw_lus[np.ix_(lat_mask, lon_mask)]
+
+
 # --------------------------------------------------------------------------
 # --------------------------------------------------------------------------
 logging.info("")
@@ -199,9 +291,11 @@ def get_fraction_per_soil_type(lu):
 logging.info("")
 logging.info("============= Calculate LU Fraction for target grid ========")
 logging.info("")
-
+start=perf_counter()
 fracs = calculate_soil_fraction(tg, soil_types, neighbor_ids, ncpu=2)
+end=perf_counter()
 
+print("Elapsed time:", end-start)
 #--------------------------------------------------------------------------
 #--------------------------------------------------------------------------
 logging.info('')
@@ -271,3 +365,4 @@ def get_fraction_per_soil_type(lu):
 logging.info('')
 logging.info('============= extpar_art_to_buffer done =======')
 logging.info('')
+

python/lib/grid_def.py

@@ -260,6 +260,9 @@ def __init__(self, gridfile):
         self.lons = np.rad2deg(self.grid.variables["clon"][:])
         self.lats = np.rad2deg(self.grid.variables["clat"][:])
 
+        self.vlons = np.rad2deg(self.grid.variables["vlon"][:])
+        self.vlats = np.rad2deg(self.grid.variables["vlat"][:])
+
     def cdo_sellonlat(self):
         '''
         create the string for the cdo option "-sellonlatbox"