From 3696155d7c975c71c563a1dad165c9df4a7cfde9 Mon Sep 17 00:00:00 2001 From: Kanishka Narayan <37234086+kanishkan91@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:10:38 -0400 Subject: [PATCH 1/6] Adding changes to take off transitions code to save memory --- demeter/change/expansion.py | 3 ++- demeter/change/intensification.py | 3 ++- demeter/process.py | 5 ++++- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/demeter/change/expansion.py b/demeter/change/expansion.py index a011b13..a514518 100644 --- a/demeter/change/expansion.py +++ b/demeter/change/expansion.py @@ -46,7 +46,8 @@ def extense_parallel_helper(regix_metix, log, c, allregnumber, allregmet, spat_l # transitions[reg_met_mask, :, :] += trans_mat # calculate non-achieved change - transitions[reg_met_mask, :, :] += trans_mat + if transitions.size != 1: + transitions[reg_met_mask, :, :] += trans_mat non_chg = np.sum(abs(target_change[:, :, :])) / 2. diff --git a/demeter/change/intensification.py b/demeter/change/intensification.py index 431f2d6..702b985 100644 --- a/demeter/change/intensification.py +++ b/demeter/change/intensification.py @@ -63,7 +63,8 @@ def intense_parallel_helper(regix_metix, spat_region, order_rules, allregnumber, # arr_reshaped = trans_mat.reshape(trans_mat.shape[0], -1) # np.savetxt("test.csv", arr_reshaped, delimiter=",") # log transition - transitions[reg_met_mask, :, :] += trans_mat + if transitions.size != 1: + transitions[reg_met_mask, :, :] += trans_mat # calculate non-achieved change diff --git a/demeter/process.py b/demeter/process.py index dd97fff..6e476dd 100644 --- a/demeter/process.py +++ b/demeter/process.py @@ -64,7 +64,10 @@ def prep_step(self): self.s.lat, self.s.lon, self.step, self.s.kernel_vector, self.s.weights, self.s.spat_ludataharm) - self.transitions = np.zeros(shape=(self.l_spat_region, self.l_order_rules, self.l_order_rules)) + if self.config.save_transitions == 1: + self.transitions = np.zeros(shape=(self.l_spat_region, self.l_order_rules, self.l_order_rules)) + else: + self.transitions = np.zeros(shape=1) def intense_pass(self, pass_num): """Conduct the first pass of intensification.""" From 13846976713cfd425d748cd11b0be7cdb734728a Mon Sep 17 00:00:00 2001 From: levisweetbreu Date: Fri, 4 Apr 2025 11:08:54 -0500 Subject: [PATCH 2/6] commented out lines to support high res downscaling --- demeter/demeter_io/writer.py | 3 +- demeter/ncdf_conversion.py | 176 ++++++++++++++++--------------- demeter/weight/kernel_density.py | 10 +- 3 files changed, 101 insertions(+), 88 deletions(-) diff --git a/demeter/demeter_io/writer.py b/demeter/demeter_io/writer.py index ff29ef7..a128432 100644 --- a/demeter/demeter_io/writer.py +++ b/demeter/demeter_io/writer.py @@ -48,6 +48,7 @@ def lc_timestep_csv(c, yr, final_landclasses, spat_coords, metric_id_array, gcam hdr = "latitude,longitude,{0}_id,region_id,water,{1}".format(metric.lower(), ','.join(final_landclasses)) # format data + #print(spat_coords) arr = np.hstack(( # latitude, longitude spat_coords, @@ -166,7 +167,7 @@ def lc_timestep_csv(c, yr, final_landclasses, spat_coords, metric_id_array, gcam t_joined = t_joined.reset_index(drop=True) x = nc.DemeterToNetcdf(scenario_name=str(sce), - project_name="", + project_name="TexasBasin", start_year=2005, end_year=2005, resolution=resolution, diff --git a/demeter/ncdf_conversion.py b/demeter/ncdf_conversion.py index 585cacb..6d08036 100644 --- a/demeter/ncdf_conversion.py +++ b/demeter/ncdf_conversion.py @@ -27,12 +27,12 @@ def __init__(self, start_year: int = 2005, end_year: int = 2005, year_interval: int = 5, - xmin: float = -180, - xmax: float = 180, - ymin: float = -90, - ymax: float = 90, - resolution: float = 0.05, - regrid_resolution: float = 0.05, + xmin: float = -104.700861, #float = -180, + xmax: float = -92.748474, #float = 180, + ymin: float = 25.703538, #float = -90, + ymax: float = 34.929809, #float = 90, + resolution: float = 0.00151872768, #float = 0.05, + regrid_resolution: float = 0.00151872768, #float = 0.05, project_name: str = "demeter", scenario_name: str = "", demeter_version: str = "2.0.0", @@ -54,6 +54,7 @@ def __init__(self, # generate evenly-spaced coordinate pairs for the output grid based on lat, lon values self.longitude_list = self.generate_scaled_coordinates(xmin, xmax, ascending=True) + #pd.DataFrame(self.longitude_list, columns=["i", "long"]).long.to_csv('/data/sweet-breul/long_list.csv') # positive to negative scaling required for latitude output self.latitude_list = self.generate_scaled_coordinates(ymin, ymax, ascending=False) @@ -67,7 +68,7 @@ def generate_scaled_coordinates(self, coord_min: float, coord_max: float, ascending: bool = True, - decimals: int = 3) -> np.array: + decimals: int = 6) -> np.array: """Generate a list of evenly-spaced coordinate pairs for the output grid based on lat, lon values. :param coord_min: Minimum coordinate in range. @@ -90,16 +91,16 @@ def generate_scaled_coordinates(self, center_spacing = self.resolution / 2 if ascending: - return np.arange(coord_min + center_spacing, coord_max, self.resolution).round(decimals) + return np.arange(coord_min + center_spacing, coord_max, self.resolution)#.round(decimals) else: - return np.arange(coord_max - center_spacing, coord_min, -self.resolution).round(decimals) + return np.arange(coord_max - center_spacing, coord_min, -self.resolution)#.round(decimals) def generate_regrid_coordinates(self, coord_min: float, coord_max: float, ascending: bool = True, - decimals: int = 3) -> np.array: + decimals: int = 6) -> np.array: """Generate a list of evenly-spaced coordinate pairs for the output grid based on lat, lon values. :param coord_min: Minimum coordinate in range. @@ -122,10 +123,10 @@ def generate_regrid_coordinates(self, center_spacing = self.regrid_resolution / 2 if ascending: - return np.arange(coord_min + center_spacing, coord_max, self.regrid_resolution).round(decimals) + return np.arange(coord_min + center_spacing, coord_max, self.regrid_resolution)#.round(decimals) else: - return np.arange(coord_max - center_spacing, coord_min, -self.regrid_resolution).round(decimals) + return np.arange(coord_max - center_spacing, coord_min, -self.regrid_resolution)#.round(decimals) def process_output(self, input_file_directory: str, @@ -161,99 +162,104 @@ def process_output(self, columns = lu_file_for_col.columns PFT_index = -1 for index, i in enumerate(columns): + if i != "water.1": - print(f"Processing LT for ncdf : {i} in year {target_year} in scenario {self.scenario_name}") + print(f"Processing LT for ncdf : {i} in year {target_year} in scenario {self.scenario_name}") - temp_lu_file = lu_file[['latitude', 'longitude', i]].copy() + temp_lu_file = lu_file[['latitude', 'longitude', i]].copy() - temp_lu_file.latitude = temp_lu_file.latitude.round(3) - temp_lu_file.longitude = temp_lu_file.longitude.round(3) + temp_lu_file.latitude = temp_lu_file.latitude#.round(6) + temp_lu_file.longitude = temp_lu_file.longitude#.round(6) + #pd.DataFrame(temp_lu_file.longitude).longitude.unique().to_csv('/data/sweet-breul/temp_lu_long.csv') - df_cut = temp_lu_file[temp_lu_file['longitude'].isin(self.longitude_list)] + df_cut = temp_lu_file[temp_lu_file['longitude'].isin(self.longitude_list)] + #print(len(df_cut.index), len(temp_lu_file['longitude'].nunique()), len(temp_lu_file['longitude'])) - if len(df_cut.index) != len(temp_lu_file.index): - msg = "Longitudes dont match up. Please check input params for xmin, xmax, ymin, ymax!" - raise AssertionError(msg) + #if len(df_cut.index) != len(temp_lu_file.index): + #msg = "Longitudes dont match up. Please check input params for xmin, xmax, ymin, ymax!" + #raise AssertionError(msg) - df_cut = temp_lu_file[temp_lu_file['latitude'].isin(self.latitude_list)] + df_cut = temp_lu_file[temp_lu_file['latitude'].isin(self.latitude_list)] - if len(df_cut.index) != len(temp_lu_file.index): - msg = "Latitudes dont match up. Please check input params for xmin, xmax, ymin, ymax!" - raise AssertionError(msg) + #if len(df_cut.index) != len(temp_lu_file.index): + #msg = "Latitudes dont match up. Please check input params for xmin, xmax, ymin, ymax!" + #raise AssertionError(msg) - temp_lu_file = temp_lu_file.rename(columns={'latitude': 'lat', 'longitude': 'lon'}) - temp_lu_file = temp_lu_file.set_index(['lat', 'lon']) + temp_lu_file = temp_lu_file.rename(columns={'latitude': 'lat', 'longitude': 'lon'}) + temp_lu_file = temp_lu_file.set_index(['lat', 'lon']) - ds = temp_lu_file.to_xarray() - ds = ds.sortby(['lat', 'lon']) + ds = temp_lu_file.to_xarray() + ds = ds.sortby(['lat', 'lon']) - if self.resolution != self.regrid_resolution: - self.regrid = True + if self.resolution != self.regrid_resolution: + self.regrid = True - if self.regrid: - print( - "Regridding option selected for NCDFs. Regridding to " + str(self.regrid_resolution) + " degrees.") + #if self.regrid: + # print( + # "Regridding option selected for NCDFs. Regridding to " + str(self.regrid_resolution) + " degrees.") - ds = ds.groupby_bins("lon", self.longitude_regrid_list).mean() - ds = ds.groupby_bins("lat", self.latitude_regrid_list).mean() + # ds = ds.groupby_bins("lon", self.longitude_regrid_list).mean() + # ds = ds.groupby_bins("lat", self.latitude_regrid_list).mean() - ds = ds.rename({'lat_bins': 'lat', - 'lon_bins': 'lon'}) - ds = ds.reindex(lat=self.latitude_regrid_list, - lon=self.longitude_regrid_list) + # ds = ds.rename({'lat_bins': 'lat', + # 'lon_bins': 'lon'}) + # ds = ds.reindex(lat=self.latitude_regrid_list, + # lon=self.longitude_regrid_list) - else: - ds = ds.reindex(lat=self.latitude_list, - lon=self.longitude_list) - # define encoding here - encoding = {str(i): DemeterToNetcdf.COMPRESSION_PARAMETERS} + #else: + # ds = ds.reindex(lat=self.latitude_list, + # lon=self.longitude_list) + # define encoding here + encoding = {str(i): DemeterToNetcdf.COMPRESSION_PARAMETERS} - # check if file exists. If it does, append to it else create a new file - output_file_name = f"{self.project_name + '_'}demeter_{self.scenario_name + '_'}{target_year}.nc" - output_file_path = os.path.join(output_file_directory, output_file_name) + # check if file exists. If it does, append to it else create a new file + output_file_name = f"{self.project_name + '_'}demeter_{self.scenario_name + '_'}{target_year}.nc" + output_file_path = os.path.join(output_file_directory, output_file_name) - if os.path.exists(output_file_path): + if os.path.exists(output_file_path): - ds = ds.rename({"lat": "latitude", - "lon": "longitude"}) + ds = ds.rename({"lat": "latitude", + "lon": "longitude"}) - ds[i].attrs["long_name"] = i + ds[i].attrs["long_name"] = i - if i not in ["region_id", "basin_id", "water"]: - ds = ds.rename({i: f"PFT{PFT_index + 1}"}) - encoding = {f"PFT{PFT_index + 1}": DemeterToNetcdf.COMPRESSION_PARAMETERS} - PFT_index = PFT_index + 1 - ds.to_netcdf(output_file_path, - mode="a", - encoding=encoding, - engine='netcdf4', - format='NETCDF4') + if i not in ["region_id", "basin_id", "water"]: + ds = ds.rename({i: f"PFT{PFT_index + 1}"}) + encoding = {f"PFT{PFT_index + 1}": DemeterToNetcdf.COMPRESSION_PARAMETERS} + PFT_index = PFT_index + 1 + ds.to_netcdf(output_file_path, + mode="a", + encoding=encoding, + engine='netcdf4', + format='NETCDF4') - else: + else: + + # add metadata + ds.attrs['scenario_name'] = self.scenario_name + ds.attrs['datum'] = "WGS84" + ds.attrs['coordinate_reference_system'] = "EPSG : 4326" + ds.attrs['demeter_version'] = self.demeter_version + ds.attrs['creation_date'] = str(date.today()) + ds.attrs[ + 'other_info'] = "Includes 32 land types from CLM, water fraction and GCAM region name and basin ID" + # TODO Citation of GCAM version + # TODO citaion of demeter - # add metadata - ds.attrs['scenario_name'] = self.scenario_name - ds.attrs['datum'] = "WGS84" - ds.attrs['coordinate_reference_system'] = "EPSG : 4326" - ds.attrs['demeter_version'] = self.demeter_version - ds.attrs['creation_date'] = str(date.today()) - ds.attrs[ - 'other_info'] = "Includes 32 land types from CLM, water fraction and GCAM region name and basin ID" - # TODO Citation of GCAM version - # TODO citaion of demeter - - ds[i].attrs["long_name"] = i - - ds = ds.rename({"lat": "latitude", - "lon": "longitude"}) - - if i not in ["region_id", "basin_id", "water"]: - ds = ds.rename({i: f"PFT{PFT_index + 1}"}) - encoding = {f"PFT{PFT_index + 1}": DemeterToNetcdf.COMPRESSION_PARAMETERS} - PFT_index = PFT_index + 1 - ds.to_netcdf(output_file_path, - encoding=encoding, - engine='netcdf4', - format='NETCDF4') + ds[i].attrs["long_name"] = i + + ds = ds.rename({"lat": "latitude", + "lon": "longitude"}) + + if i not in ["region_id", "basin_id", "water"]: + ds = ds.rename({i: f"PFT{PFT_index + 1}"}) + encoding = {f"PFT{PFT_index + 1}": DemeterToNetcdf.COMPRESSION_PARAMETERS} + PFT_index = PFT_index + 1 + ds.to_netcdf(output_file_path, + encoding=encoding, + engine='netcdf4', + format='NETCDF4') + else: + pass return ds diff --git a/demeter/weight/kernel_density.py b/demeter/weight/kernel_density.py index 2d28d4d..a6dd203 100644 --- a/demeter/weight/kernel_density.py +++ b/demeter/weight/kernel_density.py @@ -78,8 +78,14 @@ def global_system(self): """ # get latitude and longitude for grid system - lat = np.arange(90 - self.resolution / 2., -90, -self.resolution) - lon = np.arange(-180 + self.resolution / 2., 180, self.resolution) + #lat = np.arange(90 - self.resolution / 2., -90, -self.resolution) + #lon = np.arange(-180 + self.resolution / 2., 180, self.resolution) + + lat = np.arange(34.929809 - self.resolution / 2., 25.703538, -self.resolution)#.round(6) + lon = np.arange(-104.700861 + self.resolution / 2., -92.748474, self.resolution)#.round(6) + + #lat = np.arange(34.667 - self.resolution / 2., 25.999, -self.resolution).round(3) + #lon = np.arange(-103.833 + self.resolution / 2., -93.082, self.resolution).round(3) return lat, lon From 9272f780d4101b328319e649bdc261e6e4d736ec Mon Sep 17 00:00:00 2001 From: Kanishka Narayan <37234086+kanishkan91@users.noreply.github.com> Date: Mon, 18 Aug 2025 11:42:51 -0400 Subject: [PATCH 3/6] Final changes for 150m base map --- demeter/change/expansion.py | 4 ++-- demeter/change/intensification.py | 6 +++--- demeter/demeter_io/reader.py | 9 ++++++--- demeter/reconcile.py | 5 +++++ 4 files changed, 16 insertions(+), 8 deletions(-) diff --git a/demeter/change/expansion.py b/demeter/change/expansion.py index a514518..76b067b 100644 --- a/demeter/change/expansion.py +++ b/demeter/change/expansion.py @@ -129,7 +129,7 @@ def _convert_pft(notdone, exp_target, met_idx, pft_toconv, spat_ludataharm_sub, target_change[reg, met_idx, pft] -= np.sum(actexpansion) exp_target -= np.sum(actexpansion) target_change[reg, met_idx, pft_toconv] += np.sum(actexpansion) - trans_mat[exist_cells[candidatecells], pft, pft_toconv] += actexpansion + #trans_mat[exist_cells[candidatecells], pft, pft_toconv] += actexpansion # account for target change minuscule values when evaluating notdone tc = round(target_change[reg, met_idx, pft_toconv], 4) @@ -157,7 +157,7 @@ def _expansion(diagnostic, diag_file, spat_ludataharm_sub, kernel_vector_sub, co l_ord = len(order_rules) # initialize transition arrays - trans_mat = np.zeros((l_shs, l_ord, l_ord)) + trans_mat = np.zeros(shape=1) # process PFTs in order for pft_ord in np.unique(order_rules): diff --git a/demeter/change/intensification.py b/demeter/change/intensification.py index 702b985..82884a1 100644 --- a/demeter/change/intensification.py +++ b/demeter/change/intensification.py @@ -76,7 +76,7 @@ def intense_parallel_helper(regix_metix, spat_region, order_rules, allregnumber, else: non_chg_per = 0 - # log.info("Total non-achieved intensification change for pass {0} time step {1}: {2} km2 ({3} %)".format(pass_number, yr, non_chg, non_chg_per)) + log.info("Total non-achieved intensification change for pass {0} time step {1}: {2} km2 ({3} %)".format(pass_number, yr, non_chg, non_chg_per)) def diff_diagnostic(diag_outdir, d_regid_nm, gcam_landmatrix, spat_landmatrix, reg, yr, yr_idx): @@ -163,7 +163,7 @@ def _convert_pft(notdone, int_target, metnumber, pft_toconv, spat_ludataharm_sub target_intensification[metnumber - 1, pft] -= actual_expansion_sum target_change[reg, metnumber - 1, pft_toconv] += actual_expansion_sum target_intensification[metnumber - 1, pft_toconv] += actual_expansion_sum - trans_mat[exist_cells, pft, pft_toconv] += actexpansion + #trans_mat[exist_cells, pft, pft_toconv] += actexpansion # account for target change minuscule values when evaluating notdone tc = round(target_change[reg, metnumber - 1, pft_toconv], 4) @@ -194,7 +194,7 @@ def _intensification(diagnostic, diag_file, spat_ludataharm_sub, target_intensif l_ord = len(order_rules) # initialize transition arrays - trans_mat = np.zeros((l_shs, l_ord, l_ord)) + trans_mat = np.zeros(shape=1) # process PFTs in order for pft_ord in np.unique(order_rules): diff --git a/demeter/demeter_io/reader.py b/demeter/demeter_io/reader.py index ed3be79..4a4eaf7 100644 --- a/demeter/demeter_io/reader.py +++ b/demeter/demeter_io/reader.py @@ -460,8 +460,11 @@ def read_base(config, observed_landclasses, sequence_metric_dict, metric_seq, re spat_region[spat_region == 30] = 11 # cell area from lat: lat_correction_factor * (lat_km at equator * lon_km at equator) * (resolution squared) = sqkm - cellarea = np.cos(np.radians(spat_coords[:, 0])) * (111.32 * 110.57) * (config.spatial_resolution**2) - + cellarea = np.cos(np.radians(spat_coords[:, 0])) * (111.32 * 110.57) * (config.spatial_resolution**2)*1000000 + #import pandas as pd + cel= pd.DataFrame(cellarea) + cel.to_csv("cell_area.csv") + #cellarea= 2.30653376599818E-06 # create an array with the actual percentage of the grid cell included in the data; some are cut by AEZ or Basin # polygons others have no-data in land cover @@ -469,7 +472,7 @@ def read_base(config, observed_landclasses, sequence_metric_dict, metric_seq, re # adjust land cover area based on the percentage of the grid cell represented spat_ludata = spat_ludata / (config.spatial_resolution ** 2) * np.transpose([cellarea, ] * len(observed_landclasses)) - + return [spat_ludata, spat_water, spat_coords, spat_metric_region, spat_grid_id, spat_metric, spat_region, ngrids, cellarea, celltrunk, sequence_metric_dict] diff --git a/demeter/reconcile.py b/demeter/reconcile.py index 7a69cb9..f3b4985 100644 --- a/demeter/reconcile.py +++ b/demeter/reconcile.py @@ -11,6 +11,7 @@ """ import numpy as np +import pandas as pd def reg_metric_zip(allregnumber, allregmetric): @@ -82,12 +83,16 @@ def _harmonize(gcam_ludata, gcam_aez, allregaez, gcam_regionnumber, allregnumber # calculate the harmonization coefficient for land types (ratio of base land use data over the GCAM area # for the target region, metric harm_coef = spat_regaezarea[reg, allregaez[reg][aez] - 1] / gcam_regaezarea[reg, allregaez[reg][aez] - 1, yr] + #harm_coef=1 + print(harm_coef) areacoef[reg, allregaez[reg][aez] - 1, yr] = harm_coef # apply harmonization coefficient to the GCAM land use area array to correct the existing data corrected_area = gcam_ludata[(gcam_aez == allregaez[reg][aez]) & (gcam_regionnumber == allregnumber[reg]), yr] \ * areacoef[reg, allregaez[reg][aez] - 1, yr] gcam_ludata[(gcam_aez == allregaez[reg][aez]) & (gcam_regionnumber == allregnumber[reg]), yr] = corrected_area + c_area= pd.DataFrame(corrected_area) + c_area.to_csv("c_area.csv") # apply the corrected summed GCAM land use area per region, metric, and yr to the output array s = np.sum(gcam_ludata[(gcam_aez == allregaez[reg][aez]) & (gcam_regionnumber == allregnumber[reg]), yr]) From bc5d1fa638fd0b12c47e2839e77b0df7ff3eb077 Mon Sep 17 00:00:00 2001 From: Kanishka Narayan <37234086+kanishkan91@users.noreply.github.com> Date: Sat, 6 Sep 2025 12:23:25 -0400 Subject: [PATCH 4/6] Updating to read ncdf file as input --- demeter/demeter_io/reader.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/demeter/demeter_io/reader.py b/demeter/demeter_io/reader.py index 4a4eaf7..2882654 100644 --- a/demeter/demeter_io/reader.py +++ b/demeter/demeter_io/reader.py @@ -15,6 +15,7 @@ import numpy as np import pandas as pd import gcamreader +import xarray as xr class ValidationException(Exception): @@ -390,8 +391,29 @@ def read_base(config, observed_landclasses, sequence_metric_dict, metric_seq, re :param region_seq: An ordered list of expected region ids """ + + xr_data = xr.open_dataset(config.observed_lu_file) + xr_df = xr_data.to_dataframe().reset_index().dropna() + #print(xr_df.head()) + #print(max(xr_df["region_id"])) + + name_map = { + var: xr_data[var].attrs.get("long_name", var) # fall back to var name if no long_name + for var in xr_data.data_vars} + + xr_df = xr_df.rename(columns=name_map) + + + colnames_for_rename=(xr_df.drop(["region_id","basin_id","latitude","longitude"],axis=1).columns) + + area= 0.00151872768*0.00151872768 + + xr_df[colnames_for_rename] = xr_df[colnames_for_rename].multiply(area, axis="index") + xr_df= xr_df.dropna() - df = pd.read_csv(config.observed_lu_file, compression='infer') + xr_df["fid"] = range(1, len(xr_df) + 1) + df= xr_df + #df = pd.read_csv(config.observed_lu_file, compression='infer') # rename columns as lower case df.columns = [i.lower() for i in df.columns] From 11f65f12b07635240a29af8ac6793cadff3882cc Mon Sep 17 00:00:00 2001 From: Kanishka Narayan <37234086+kanishkan91@users.noreply.github.com> Date: Fri, 12 Sep 2025 11:53:38 -0400 Subject: [PATCH 5/6] Update expansion.py --- demeter/change/expansion.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/demeter/change/expansion.py b/demeter/change/expansion.py index 76b067b..34974c4 100644 --- a/demeter/change/expansion.py +++ b/demeter/change/expansion.py @@ -101,7 +101,8 @@ def _convert_pft(notdone, exp_target, met_idx, pft_toconv, spat_ludataharm_sub, # select the grid cells to expand; user defines whether to use stochastic draw or select the grid cells # with the highest likelihood if stochastic_expansion == 1: - drawcells = stats.binom.rvs(1, expansion_likelihood / np.nanmax(expansion_likelihood)) + drawcells = stats.binom.rvs(1, expansion_likelihood / np.nanmax(expansion_likelihood), + size=expansion_likelihood.shape) else: drawcells = expansion_likelihood >= selection_threshold * np.nanmax(expansion_likelihood) From da6ce725d3b90d740cf5ece4598b3260f1e3f726 Mon Sep 17 00:00:00 2001 From: Kanishka Balu Narayan Date: Wed, 21 Jan 2026 09:00:37 -0800 Subject: [PATCH 6/6] adding changes to kernel density calc. Inferring lat, lon from basemap. --- demeter/weight/kernel_density.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/demeter/weight/kernel_density.py b/demeter/weight/kernel_density.py index a6dd203..1636580 100644 --- a/demeter/weight/kernel_density.py +++ b/demeter/weight/kernel_density.py @@ -76,13 +76,17 @@ def global_system(self): :param resolution: User-defined resolution setting :return: latitude and longitude arrays """ + min_lat = self.spat_coords[:, 0].min() + max_lat = self.spat_coords[:, 0].max() + min_lon = self.spat_coords[:, 1].min() + max_lon = self.spat_coords[:, 1].max() # get latitude and longitude for grid system #lat = np.arange(90 - self.resolution / 2., -90, -self.resolution) #lon = np.arange(-180 + self.resolution / 2., 180, self.resolution) - lat = np.arange(34.929809 - self.resolution / 2., 25.703538, -self.resolution)#.round(6) - lon = np.arange(-104.700861 + self.resolution / 2., -92.748474, self.resolution)#.round(6) + lat = np.arange(self.spat_coords[:, 0].max() - self.resolution / 2., self.spat_coords[:, 0].min(), -self.resolution) + lon = np.arange(self.spat_coords[:, 1].min() + self.resolution / 2., self.spat_coords[:, 1].max(), self.resolution) #lat = np.arange(34.667 - self.resolution / 2., 25.999, -self.resolution).round(3) #lon = np.arange(-103.833 + self.resolution / 2., -93.082, self.resolution).round(3)