Script 2 — Baseline calculation
The core of the procedure: the spectral unmixing and the stocking index, the screening, the matching of control to project plots, the match quality, the slopes and the ex-post performance benchmark.
The listing below is vm_0047_baseline_calculation_v2.py: 708 lines, the script delivered earlier with the changes of Table 19, each marked v2 in the code. It is read from the file at render time, so the document and the script cannot drift apart.
⬇ Download vm_0047_baseline_calculation_v2.py
#(c) Terra Global Capital CONFIDENTIAL
import geopandas as gpd
import pandas as pd
import glob
import pandas as pd
import numpy as np
from itertools import combinations
import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF
from scipy.optimize import linear_sum_assignment
from scipy.spatial.distance import cdist
from scipy.stats import linregress
import itertools
from itertools import permutations
import numpy as np
import pandas as pd
from sklearn.decomposition import NMF
from scipy.optimize import nnls
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from numpy.linalg import lstsq
from typing import Tuple
import os # v2: paths relative to the update folder
import rasterio # v2: MapBiomas 2018 land-cover screen
#from munkres import Munkres
# v2: NDFI of Souza et al. (2005) with the five fixed CODED endmembers (Bullock et al. 2020), in place of the NMF
# abundances. Rows: GV, NPV, Soil, Shade, Cloud; columns: Landsat 8 SR_B2-SR_B7 (B0-B5), reflectance.
CODED_ENDMEMBERS = np.array([[0.05, 0.09, 0.04, 0.61, 0.30, 0.10],
[0.14, 0.17, 0.22, 0.30, 0.55, 0.30],
[0.20, 0.30, 0.34, 0.58, 0.60, 0.58],
[0.00, 0.00, 0.00, 0.00, 0.00, 0.00],
[0.90, 0.96, 0.80, 0.78, 0.72, 0.65]])
# v2: USGS Collection 2 Level 2 surface reflectance valid range (DN)
VALID_DN = (7273, 43636)
def souza_ndfi(dn):
# v2: fully constrained unmixing (non-negative, sum to one through a weighted row of ones) of the reflectance
# DN x 2.75e-5 - 0.2; NDFI = (GVs - (NPV + Soil)) / (GVs + NPV + Soil), GVs = GV / (1 - Shade); cloud stays in
# its own fraction and does not enter the NDFI
A = np.vstack([CODED_ENDMEMBERS.T, 1e3 * np.ones(len(CODED_ENDMEMBERS))])
refl = np.asarray(dn, dtype=float) * 2.75e-5 - 0.2
F = np.array([nnls(A, np.append(r, 1e3))[0] for r in refl])
gv, npv, soil, shade = F[:, 0], F[:, 1], F[:, 2], F[:, 3]
with np.errstate(divide='ignore', invalid='ignore'):
gvs = gv / (1 - shade)
return (gvs - (npv + soil)) / (gvs + npv + soil)
def nmf_spectral_unmixing(df_, n_components=3, custom_endmembers=None):
"""
Apply NMF (Non-negative Matrix Factorization) for spectral unmixing with an option to use custom endmembers.
:param df_: DataFrame with columns for different spectral bands.
:param n_components: Number of components (endmembers) to extract.
:param custom_endmembers: Optional. DataFrame with custom endmembers. Must have `n_components` columns.
:return: DataFrame with additional columns for the fractions of each component, and the endmember spectra.
"""
if custom_endmembers is not None and custom_endmembers.shape[0] != n_components:
raise ValueError("custom_endmembers must have the same number of columns as n_components")
if custom_endmembers is not None:
# Use custom endmembers for prediction
model = NMF(n_components=n_components, init='random', random_state=0, max_iter=1000)
H = custom_endmembers # Transpose to match the shape expected for multiplication
W_, _, _, _ = np.linalg.lstsq(H.T, df_.values.T, rcond=None) # Solve for W in H*W = df_.T
W= model.fit_transform(df_, W=W_.T, H=H.T)
else:
# Fit Non-Negative Matrix Factorization (NMF) to find endmembers
model = NMF(n_components=n_components, init='random', random_state=0, max_iter=1000)
W = model.fit_transform(df_) # Abundance fractions
H = model.components_ # Endmember spectra
# Add fractions back to the DataFrame
for i in range(n_components):
df_[f'band_{i+1}_nmf'] = W[:, i]
return df_, H
# Calculate the weighted average, excluding NaN values
def weighted_average(row, ndfi_cols, dist_cols):
ndfi_values = row[ndfi_cols].values
weights = row[dist_cols].values
# Filter out NaN values and their corresponding weights
mask = ~np.isnan(ndfi_values)
filtered_ndfi_values = ndfi_values[mask]
filtered_weights = weights[mask]
# If all values are NaN, return NaN
if len(filtered_ndfi_values) == 0:
return np.nan
# Calculate weighted average
return np.average(filtered_ndfi_values, weights=filtered_weights)
def calculate_ndfi(row):
return (row['band_1_nmf'] - (row['band_2_nmf'] + row['band_3_nmf'])) / (row['band_1_nmf'] + row['band_2_nmf'] + row['band_3_nmf'])
def get_endmembers(df_agb, n_endmenbers=3):
#df_agb = df_agb[[col for col in df_agb.columns if 'plotId' in col or 'DSM' in col or 'landcover' in col or 'B' in col or 'plot_type' in col or 'Year' in col]]
# define features as all columns but plot_type or Year
landsat_bands = [col for col in df_agb.columns if 'B' in col and col != 'AGB' ]
# remove rows where sum of columsn with band is 0
df_agb = df_agb.loc[(df_agb[landsat_bands] != 0).all(axis=1)]
df_agb[landsat_bands] = df_agb[landsat_bands].multiply(2.75e-5)
# apply l2 norm to the bands
df_agb[landsat_bands] = df_agb[landsat_bands].apply(lambda x: x / np.linalg.norm(x), axis=1)
# count nan in each column
print(df.isna().sum())
# remove rows with nan
df_agb = df_agb.dropna()
# reset index
df_agb = df_agb.reset_index(drop=True)
print(df_agb.isna().sum())
# Use NMF to get the endmembers
data_for_ndfi, endmember_spectra = nmf_spectral_unmixing(df_agb[landsat_bands], n_components=n_endmenbers, custom_endmembers=None)
# calculate NDFI from all combinations of endmembers and use the one that best correlates with the agb
# List of all column names you want to consider for combinations, exclude AGB from this list
columns_for_combination = [col for col in data_for_ndfi.columns if 'nmf' in col]
# Calculate NDFI for each combination and find the one with the highest correlation with AGB
best_combination = None
highest_correlation = -np.inf # Start with the lowest possible correlation
for combo in permutations(columns_for_combination, 3):
print(combo)
tmp_df = data_for_ndfi[list(combo)]
NDFI = (tmp_df.iloc[:, 0] - (tmp_df.iloc[:, 1] + tmp_df.iloc[:, 2])) / \
(tmp_df.iloc[:, 0] + tmp_df.iloc[:, 1] + tmp_df.iloc[:, 2])
# Calculate Pearson correlation coefficient with AGB
correlation_with_AGB = NDFI.corr(df_agb['AGB'])
print(correlation_with_AGB)
# Update best combination if this one is better
if correlation_with_AGB > highest_correlation:
highest_correlation = correlation_with_AGB
best_combination = combo
# Output the best combination and its correlation with AGB
print(fr"best combination is {best_combination}, with a correlation with AGB of: {highest_correlation}")
# calculate the NDFI using the best combination
return endmember_spectra, best_combination, highest_correlation
# Function to fill the distance column based on the outside_index_type
# Extract the index number from 'outside_index_type' and select the correct distance column dynamically
def get_dynamic_distance(row):
# Extract the number at the end of 'outside_index_type'
index_number = int(row['outside_index_type'].split('_')[-1])
# Construct the column name for the corresponding distance
dist_column = f'outside_dist_{index_number}'
# Return the value from the appropriate column
return row[dist_column]
def analyze(df, historic_threshold, endmember_spectra, best_combination, n_endmenbers=5, n_controls_per_plot = 3, out_path = '.', filename = 'summary', get_the_smallest_distance = 45, gpkg_dir = None, source_gpkg = None):
# Main analysis function to orchestrate the workflow
# TODO make this more generalizable and make a function
df = df[[col for col in df.columns if 'plotId' in col or 'DSM' in col or 'landcover' in col or 'B' in col or 'plot_type' in col or 'Year' in col]]
# define features as all columns but plot_type or Year
features = [col for col in df.columns if 'DSM' in col or 'landcover' in col or 'B' ]
landsat_bands = [col for col in df.columns if 'B' in col ]
# remove rows where sum of columsn with band is 0
df = df.loc[(df[features] != 0).all(axis=1)]
# v2: keep the DN for the NDFI and the valid-range mask; reflectance with the -0.2 offset before the L2 norm
dn = df[landsat_bands].astype(float)
df[landsat_bands] = df[landsat_bands].multiply(2.75e-5) - 0.2
# apply l2 norm to the bands
df[landsat_bands] = df[landsat_bands].apply(lambda x: x / np.linalg.norm(x), axis=1)
# count nan in each column
print(df.isna().sum())
# remove rows with nan
df = df.dropna()
# v2: Souza NDFI from the DN, in place of the NMF; a plot with any row outside the valid DN range or with an
# undefined NDFI is dropped (the code's rule for missing data)
df['band_NDFI'] = souza_ndfi(dn.loc[df.index].values)
valid_dn = ((dn.loc[df.index] >= VALID_DN[0]) & (dn.loc[df.index] <= VALID_DN[1])).all(axis=1).values
bad_plots = df.loc[~valid_dn | ~np.isfinite(df['band_NDFI'].values), 'plotId'].unique()
print(f'v2: rows outside the valid DN range {(~valid_dn).sum()}, undefined NDFI {(~np.isfinite(df.band_NDFI)).sum()}, plots dropped {len(bad_plots)}')
df = df[~df.plotId.isin(bad_plots)]
df = df.reset_index(drop=True)
# Get unique plotId values and create a mapping dictionary
unique_plot_ids = df['plotId'].unique()
plot_id_mapping = {old_id: new_id for new_id, old_id in enumerate(unique_plot_ids)}
# Remap plotId column using the mapping dictionary
df['plotId'] = df['plotId'].map(plot_id_mapping)
# use the same mapping for the plotId in the source_gpkg
source_gpkg['plotId'] = source_gpkg['plotId'].map(plot_id_mapping)
df_historic = df[df.Year <= historic_threshold]
df_monitoring = df[df.Year > historic_threshold]
# Pivot the DataFrame into wide format using 'plotid' as the index
wide_df_historic = df_historic.pivot(index='plotId', columns='Year')
wide_df_monitoring = df_monitoring.pivot(index='plotId', columns='Year')
# Reset the index if you want to turn the multi-level index into columns
wide_df_historic.reset_index(inplace=True)
wide_df_monitoring.reset_index(inplace=True)
# count nan in each column
print(wide_df_historic.isna().sum())
# Flatten the column MultiIndex
wide_df_historic.columns = ['_'.join(str(col) for col in cols_tuple) if cols_tuple[1] != '' else cols_tuple[0] for cols_tuple in wide_df_historic.columns.values]
wide_df_monitoring.columns = ['_'.join(str(col) for col in cols_tuple) if cols_tuple[1] != '' else cols_tuple[0] for cols_tuple in wide_df_monitoring.columns.values]
# Separate the data
within_df_historic = wide_df_historic[wide_df_historic[f'plot_type_{ df.Year.min()}'] == 'within']
outside_df_historic = wide_df_historic[wide_df_historic[f'plot_type_{ df.Year.min()}'] == 'outside']
# Separate the data
wide_df_monitoring = wide_df_monitoring[wide_df_monitoring[f'plot_type_{df.Year.max()}'] == 'within']
wide_df_monitoring = wide_df_monitoring[wide_df_monitoring[f'plot_type_{df.Year.max()}'] == 'outside']
# Initialize an empty DataFrame to store pairs
#bands are all the columns that contain "band" in the name
bands = [col for col in wide_df_historic.columns if 'B' in col or 'DSM_2015' in col or 'topoDiversity_constant_2014' in col or 'band' in col or 'GV' in col or 'NPV' in col or 'Soil' in col]
# v2: the project plots' within_index (row position) is used below as their plotId; it must be
assert (within_df_historic['plotId'].values == np.arange(len(within_df_historic))).all()
# v2: covariates = SI of the historic years + principal components of the six 2018 bands (80 % of the
# variance), standardised over project and donor plots together, in place of the 49 covariates
n_within = len(within_df_historic)
band_cols_2018 = [col for col in bands if 'NDFI' not in col and col.endswith('_2018')]
si_cols = [col for col in bands if 'NDFI' in col]
all_plots = pd.concat([within_df_historic, outside_df_historic], ignore_index=True)
z = (all_plots[band_cols_2018] - all_plots[band_cols_2018].mean()) / all_plots[band_cols_2018].std()
_, s, vt = np.linalg.svd(z.values, full_matrices=False)
n_pcs = int(np.searchsorted(np.cumsum(s ** 2) / (s ** 2).sum(), 0.80)) + 1
pcs = [f'PC{i + 1}' for i in range(n_pcs)]
all_plots[pcs] = z.values @ vt[:n_pcs].T
bands = si_cols + pcs
all_plots[bands] = (all_plots[bands] - all_plots[bands].mean()) / all_plots[bands].std()
within_df_historic = all_plots.iloc[:n_within].reset_index(drop=True)
outside_df_historic = all_plots.iloc[n_within:].reset_index(drop=True)
# run hungarian algorithm using manhattan distance and linear sum assignment
# Calculate the initial distance matrix
distance_matrix = cdist(within_df_historic[bands], outside_df_historic[bands], metric='cityblock')
# Placeholder for results
assignments = []
# Perform iterative assignment
for _ in range(4): # assign each 'within' row to 4 'outside' rows
row_ind, col_ind = linear_sum_assignment(distance_matrix)
# Store assignments and remove these 'outside' rows from consideration
for r, c in zip(row_ind, col_ind):
assignments.append((r, c)) # Track which 'within' row is assigned to which 'outside' row
distance_matrix[:, c] = np.inf # Set distances to infinity for assigned 'outside' rows to exclude them from future iterations
distance_matrix = cdist(within_df_historic[bands], outside_df_historic[bands], metric='cityblock')
assignment_dict = []
for r, c in assignments:
distance = distance_matrix[r, c]
# Capture the distance for the current assignment
assignment_dict.append((r, c, distance))
assignment_dict_with_distances = {i: [] for i in range(len(within_df_historic))}
for r, c, distance in assignment_dict:
assignment_dict_with_distances[r].append((c, distance))
# Creating a list of lists for the DataFrame, now including distances
data_for_df_with_distances = []
for within_index in sorted(assignment_dict_with_distances.keys()):
row_data = [within_index] # Start with the 'within' index
for outside_index, distance in assignment_dict_with_distances[within_index]:
row_data.extend([outside_index, distance]) # Add 'outside' index and corresponding distance
data_for_df_with_distances.append(row_data)
# Define column names dynamically based on the number of assignments (e.g., 4 matches)
column_names = ['within_index']
for i in range(1, 5): # Assuming 4 matches
column_names.extend([f'outside_index_{i}', f'outside_dist_{i}'])
# Convert to DataFrame
pairs_df = pd.DataFrame(data_for_df_with_distances, columns=column_names)
ndfi_per_plot_per_Year = df[['plotId', 'Year', 'band_NDFI']]
if get_the_smallest_distance is not None:
# Select the 45 'within' plots with the smallest average distance to 'outside' plots
# 1. Group by 'within_index' and calculate the average 'outside_dist'
# get a list of elements in the columns oustide_index_1 to outside_index_4
long_df = pd.wide_to_long(pairs_df, stubnames=['outside_index', 'outside_dist'], i='within_index', j='number', sep='_')
long_df = long_df.reset_index()
# Drop the 'number' column since it's no longer needed
long_df = long_df.drop(columns=['number'])
print(long_df)
long_df.within_index.unique()
grouped_df = long_df.groupby('within_index')['outside_dist'].mean().reset_index()
# 2. Sort by the average 'outside_dist' in ascending order
sorted_df = grouped_df.sort_values(by='outside_dist')
# 3. Select the top 45 'within_index' with the lowest averages
top_45_within_index = sorted_df.head(get_the_smallest_distance)['within_index']
print(top_45_within_index)
# from long df, get only the rows where within_index is in top_45_within_index
pairs_df_top_45 = long_df[long_df['within_index'].isin(top_45_within_index)]
# create a list of unique within_index and outside_index from pairs_df_top_45
unique_within_index = pairs_df_top_45['within_index'].unique()
unique_outside_index = pairs_df_top_45['outside_index'].unique()
donor_df_reindexed = outside_df_historic
donor_df_reindexed=donor_df_reindexed.reset_index()
donor_df_reindexed['outside_index'] = donor_df_reindexed.index
stele_donors = donor_df_reindexed[['plotId', 'outside_index']][donor_df_reindexed.index.isin( pairs_df_top_45['outside_index'])]
uniques_wihtin_out = np.unique(np.concatenate((unique_within_index, stele_donors['plotId']), axis=None))
dataset = source_gpkg
dataset["PLOT_GROUP"] = "Donor"
dataset["DOWNSTREAM"] = False
dataset.loc[dataset["plotId"].isin(uniques_wihtin_out), "DOWNSTREAM"] = True
dataset.loc[dataset.plot_type == 'within', 'PLOT_GROUP'] = 'Project'
# create a sample of the data with only data from one year
sample = dataset[dataset.Year == 2018]
#save all_data and sample to 2 different geopavkages
dataset.to_file(fr'{gpkg_dir}\all_data_new2.gpkg', driver='GPKG')
sample.to_file(fr'{gpkg_dir}\sample_new2.gpkg', driver='GPKG')
# order pairs_df_top_45 by within_index and outside_dist
pairs_df_top_45 = pairs_df_top_45.sort_values(by=['within_index', 'outside_dist'])
# paste the plotId of the outside_index to the pairs_df_top_45
# use stele donor to join plotId to outside_index
pairs_df_top_45 = pd.merge(pairs_df_top_45, stele_donors, on='outside_index')
# save unique_within_index and unique_outside_index to csv
pd.DataFrame(pairs_df_top_45).to_csv(fr'{gpkg_dir}\lookup_and_distances_new.csv')
# keep the pairs_df within_index that are in top_45_within_index
pairs_df = pairs_df[pairs_df.within_index.isin(top_45_within_index)]
# Reshape pairs_df for easier merging
index_cols = ['within_index']
dist_cols = [col for col in pairs_df.columns if 'outside_dist' in col]
index_dist_cols = index_cols + dist_cols
melted_pairs_df = pairs_df.melt(id_vars=index_dist_cols,
value_vars=[col for col in pairs_df.columns if 'outside_index' in col],
var_name='outside_index_type',
value_name='plotId')
# v2: the melted column holds the donor's row position in outside_df_historic; map it to the donor's plotId
melted_pairs_df['plotId'] = outside_df_historic['plotId'].values[melted_pairs_df['plotId'].astype(int).values]
# calculate teh SDM from within_df_historic[bands], outside_df_historic[bands]
weights = melted_pairs_df.copy()
outside_dist_cols = [col for col in melted_pairs_df.columns if 'outside_dist' in col]
weights[outside_dist_cols] = np.exp(-weights[outside_dist_cols]).div(np.exp(-weights[outside_dist_cols]).sum(axis=1), axis=0)
# Map outside_index_type to the corresponding outside_dist column
outside_dist_map = {
'outside_index_1': 'outside_dist_1',
'outside_index_2': 'outside_dist_2',
'outside_index_3': 'outside_dist_3',
'outside_index_4': 'outside_dist_4'
}
# Create the plot_weight column by applying the mapping
weights['plot_weight'] = weights.apply(lambda row: row[outside_dist_map[row['outside_index_type']]], axis=1)
# Select the desired columns for the new DataFrame
weights = weights[['within_index', 'plotId', 'plot_weight']]
# Display the new DataFrame
print(weights)
new_list = ['plotId', ', '.join(bands)]
tmp_merger = pd.merge(weights, outside_df_historic, on='plotId')
# Calculate the weighted average for each column
weighted_averages_control = {}
weighted_std_control = {}
weighted_averages_project = {}
weighted_std_project = {}
SDM = {}
# v2: Equation (A2) on the Step 3 population: the selected project plots, not all of them
selected_project = within_df_historic.iloc[sorted(pairs_df['within_index'].unique())]
for column in bands:
weighted_sum = (tmp_merger[column] * tmp_merger['plot_weight']).sum()
total_weight = tmp_merger['plot_weight'].sum()
weighted_averages_control[column] = weighted_sum / total_weight
# standard deviation
# v2: weighted by the plot weights
weighted_std_control[column] = np.sqrt((tmp_merger['plot_weight'] * (tmp_merger[column] - weighted_averages_control[column])**2).sum() / total_weight)
weighted_averages_project[column] = selected_project[column].mean()
weighted_std_project[column] = selected_project[column].std()
# calculate the SDM for the column
# v2: pooled standard deviation, sqrt((s_p^2 + s_c^2) / 2), in place of sqrt(s_p + s_c)
SDM[column] = abs(weighted_averages_project[column] - weighted_averages_control[column]) / np.sqrt((weighted_std_project[column]**2 + weighted_std_control[column]**2) / 2)
# Convert the results to a DataFrame
SDM_df = pd.DataFrame(SDM, index=['Covariate_SDM']).T
# save the SDM_df to a csv
SDM_df.to_csv(fr'{gpkg_dir}\SDM_df.csv')
# Merge with ndfi_df
merged_df = pd.merge(melted_pairs_df, ndfi_per_plot_per_Year, on='plotId')
# count nan in each column
print(merged_df.isna().sum())
# Pivot to get back to the original shape with band_NDFI and distance values
result_df = merged_df.pivot_table(index=index_dist_cols + ['Year'],
columns='outside_index_type',
values='band_NDFI').reset_index()
# Join with ndfi_df to add NDFI_within
final_df = result_df.merge(ndfi_per_plot_per_Year.rename(columns={'plotId': 'within_index', 'band_NDFI': 'NDFI_within'}), on=['within_index', 'Year'])
# rescale disntances to sum to 1
# v2: weights of Equation (2), exp(-MD) / sum(exp(-MD)), in place of the distances rescaled to sum to 1
final_df[outside_dist_cols] = np.exp(-final_df[outside_dist_cols]).div(np.exp(-final_df[outside_dist_cols]).sum(axis=1), axis=0)
print(final_df)
#save the final_df to a csv
final_df.to_csv(fr'{out_path}\weights_and_si_df.csv')
# Columns for band_NDFI values and their corresponding distances
ndfi_cols = ['outside_index_1', 'outside_index_2', 'outside_index_3', 'outside_index_4']
dist_cols = ['outside_dist_1', 'outside_dist_2', 'outside_dist_3', 'outside_dist_4']
# Apply the function to each row
final_df['weighted_avg_NDFI'] = final_df.apply(weighted_average, args=(ndfi_cols, dist_cols), axis=1)
final_df.to_csv(fr'{out_path}\final_df.csv') #from here using the last dataset of pair points
data_for_regression = final_df[['within_index', 'Year','NDFI_within', 'weighted_avg_NDFI']]
# use only the data with Year > than historic_threshold
data_for_benchmark = data_for_regression[data_for_regression.Year >= historic_threshold]
# Assuming 'data_for_benchmark' is a pandas DataFrame with 'Year', 'NDFI_within', and 'weighted_avg_NDFI' columns
# Initialize an empty DataFrame to hold the results
columns = ['t', 'SI_project', 'dSI_WPT', 'SI_control', 'dSI_Control', 'Z_Slopes', 'p_value_control', 'PBt', 'PBratio']
final_table = pd.DataFrame(columns=columns)
data_for_baseline = data_for_regression[data_for_regression.Year <= historic_threshold]
delta_si_control_history, _, _, p_value_control_history, std_err_control_history = linregress(data_for_baseline['Year'], data_for_baseline['weighted_avg_NDFI'])
# get the data from the first Year after the historic_threshold
for Year in range(historic_threshold, data_for_benchmark['Year'].max() + 1):
if Year == historic_threshold:
continue
data_for_Year = data_for_benchmark[data_for_benchmark['Year'] <= Year]
# Calculate the slope for NDFI_within and weighted_avg_NDFI
delta_si_wpt, _, _, p_value_sample, std_err_sample = linregress(data_for_Year['Year'], data_for_Year['NDFI_within'])
delta_si_control, _, _, p_value_control, std_err_control = linregress(data_for_Year['Year'], data_for_Year['weighted_avg_NDFI'])
# Calculate the significance of the difference between slopes
z_slopes = (delta_si_wpt - delta_si_control) / np.sqrt(std_err_sample**2 + std_err_control**2)
# Calculate PBt
PBt = 1
if abs(z_slopes) > 1.96:
if p_value_control > 0.05 or delta_si_control < 0:
PBt = 0
else:
PBt = min(1.0, delta_si_control / delta_si_wpt)
if PBt < 0 and delta_si_control > 0:
PBt = 1
# Append the results for this Year to the DataFrame
data_for_Year = data_for_benchmark[data_for_benchmark['Year'] == Year]
final_table = pd.concat([final_table, pd.DataFrame([{
't': Year- historic_threshold,
'SI_project': data_for_Year['NDFI_within'].mean(), # Assuming you want the average NDFI_within for the Year
'dSI_WPT': delta_si_wpt,
'SI_control': data_for_Year['weighted_avg_NDFI'].mean(), # Assuming you want the average weighted_avg_NDFI for the Year
'dSI_Control': delta_si_control,
'std_err_project': std_err_sample,
'std_err_control': std_err_control,
'Z_Slopes': z_slopes,
'p_value_control': p_value_control,
'PBt': PBt,
'PBratio': delta_si_control / delta_si_wpt # check the actual value of the ratio
}])], ignore_index=True)
# get working directoru
#final_table.to_csv(fr'{filename}.csv') # Replace with the path where you want to
print(final_table)
final_table.to_csv(fr'{out_path}\final_table.csv')
'''
print(final_table)
# plot the data for regression using Year as the x, NDFI_within and weighted_avg_NDFI as y (blue and red respectively)
# Plotting
plt.figure(figsize=(10, 6))
# Plot NDFI_within
plt.scatter(final_df['Year'], final_df['NDFI_within'], color='blue', label='NDFI_within')
# Plot weighted_avg_NDFI
plt.scatter(final_df['Year'], final_df['weighted_avg_NDFI'], color='red', label='weighted_avg_NDFI')
# Adding titles and labels
plt.title('NDFI_within and Weighted Average NDFI Over Years')
plt.xlabel('Year')
plt.ylabel('NDFI Values')
plt.legend()
plt.show()
'''
# repeat using the average across plots
#select data for benchmark in 2020
mean_benchmark = data_for_benchmark.groupby(['Year']).mean().reset_index()
# only get Years after threshold
mean_benchmark = mean_benchmark[mean_benchmark.Year >= historic_threshold]
mean_benchmark_baseline = data_for_baseline.groupby(['Year']).mean().reset_index()
delta_si_control_mean_history, _, _, p_value_control_mean_history, std_err_control_mean_history = linregress(mean_benchmark_baseline['Year'], mean_benchmark_baseline['weighted_avg_NDFI'])
final_table_avg = pd.DataFrame(columns=columns)
for Year in range(historic_threshold+1, mean_benchmark['Year'].max() +1):
data_for_Year = mean_benchmark[mean_benchmark['Year'] <= Year]
# Calculate the slope for NDFI_within and weighted_avg_NDFI
delta_si_wpt, _, _, p_value_sample, std_err_sample = linregress(data_for_Year['Year'], data_for_Year['NDFI_within'])
delta_si_control, _, _, p_value_control, std_err_control = linregress(data_for_Year['Year'], data_for_Year['weighted_avg_NDFI'])
# Calculate the significance of the difference between slopes
z_slopes = (delta_si_wpt - delta_si_control) / np.sqrt(std_err_sample**2 + std_err_control**2)
# Calculate PBt
PBt = 1
if abs(z_slopes) > 1.96:
if p_value_control > 0.05 or delta_si_control < 0:
PBt = 0
else:
PBt = min(1.0, delta_si_control / delta_si_wpt)
# Append the results for this Year to the DataFrame
#data_for_Year = mean_benchmark[mean_benchmark['Year'] == 2020]
final_table_avg = pd.concat([final_table_avg, pd.DataFrame([{
't': Year- historic_threshold,
'SI_project': data_for_Year['NDFI_within'].mean(), # Assuming you want the average NDFI_within for the Year
'dSI_WPT': delta_si_wpt,
'SI_control': data_for_Year['weighted_avg_NDFI'].mean(), # Assuming you want the average weighted_avg_NDFI for the Year
'dSI_Control': delta_si_control,
'std_err_project': std_err_sample,
'std_err_control':std_err_control,
'Z_Slopes': z_slopes,
'p_value_control': p_value_control,
'PBt': PBt,
'PBratio': delta_si_control / delta_si_wpt # Placeholder for 'PBratio
}])], ignore_index=True)
print(final_table_avg)
#final_table_avg.to_csv(fr'{filename}_avg.csv') # Replace with the path where you want to
#mean_benchmark.to_csv(fr'SI_data_average_{Year}_final.csv')
#data_for_regression.to_csv(fr'SI_data_all_points_{Year}_final.csv')
# group by Year and plotId and get the mean of the NDFI_within and weighted_avg_NDFI
final_df_grouped = final_df.groupby(['Year']).mean().reset_index()
# remove Years before 2009
#final_df_grouped = final_df_grouped[final_df_grouped.Year >= 2009]
plt.figure(figsize=(12, 6))
# Plot NDFI_within
plt.scatter(final_df_grouped['Year'], final_df_grouped['NDFI_within'], color='blue', label='SI project')
# Plot weighted_avg_NDFI
plt.scatter(final_df_grouped['Year'], final_df_grouped['weighted_avg_NDFI'], color='red', label='SI control')
# Filter the data for trendlines from threshold Year
filtered_df = final_df_grouped[final_df_grouped['Year'] >= historic_threshold+1]
# Add trendlines only from threshold Year
z_within = np.polyfit(filtered_df['Year'], filtered_df['NDFI_within'], 1)
p_within = np.poly1d(z_within)
plt.plot(filtered_df['Year'], p_within(filtered_df['Year']), "b--")
z_weighted = np.polyfit(filtered_df['Year'], filtered_df['weighted_avg_NDFI'], 1)
p_weighted = np.poly1d(z_weighted)
plt.plot(filtered_df['Year'], p_weighted(filtered_df['Year']), "r--")
# Add a vertical line at the historic_threshold
plt.axvline(x=historic_threshold+1, color='k', linestyle='--')
# Adding titles and labels
plt.title('SI project (blue) and SI control (red)')
plt.xlabel('Year')
plt.ylabel('SI (NDFI)')
plt.legend( ['Mean SI project', 'Mean SI control'])
plt.show()
#plt.savefig(fr'{filename}NDFI_timeseries_final.png')
return final_table, final_table_avg, delta_si_control_history, p_value_control_history, std_err_control_history, delta_si_control_mean_history,p_value_control_mean_history, std_err_control_mean_history
from scipy.spatial import distance_matrix
def select_project_plots(df, n_plots=50, numerical_cols=None):
# Select 30 project plots randomly
# Assuming df is your original DataFrame
df_within = df[df['plot_type'] == 'within']
df_outside = df[df['plot_type'] == 'outside']
# Select the numerical columns
if numerical_cols is None:
numerical_cols = df.columns[df.dtypes != 'object'].tolist()
# Compute distance matrix
dist_matrix = distance_matrix(df_within[numerical_cols].values, df_outside[numerical_cols].values)
# Find the minimum distance to any 'outside' plot for each 'within' plot
min_distances = np.min(dist_matrix, axis=1)
df_within['min_distance'] = min_distances
# Now aggregate these distances by 'plotID' across Years, if necessary
# Assuming 'plotID' and 'Year' as the indices you want to group by
agg_distances = df_within.groupby('plotId')['min_distance'].mean()
# Select the 45 plotIDs with the lowest average distance
selected_plotIDs = agg_distances.nsmallest(n_plots).index.tolist()
return df[df['plotId'].isin(selected_plotIDs)]
def compile_timeseries(directory_path):
# Use glob to find all .gpkg files in the directory
gpkg_files = glob.glob(directory_path)
# Initialize an empty DataFrame for appending all contents
all_data = pd.DataFrame()
# Loop through the list of GPKG files
for file_path in gpkg_files:
gdf = gpd.read_file(file_path)
all_data = pd.concat([all_data, gdf], ignore_index=True)
# Now, all_data contains the combined content from all the GPKG files
columns_to_merge = {}
for col in all_data.columns:
# Consider the last 4 letters of each column
key = col[-4:]
if key in columns_to_merge:
columns_to_merge[key].append(col)
else:
columns_to_merge[key] = [col]
# Step 2: Merge columns and calculate the average
for key, cols in columns_to_merge.items():
if len(cols) > 1: # Only merge if there are at least 2 columns with the same last 4 letters
# Calculate the mean using np.nanmean, axis=1 for row-wise operation
all_data[key + '_avg'] = np.nanmean(all_data[cols], axis=1)
# Remove the original columns
all_data.drop(cols, axis=1, inplace=True)
# based on the geometry column, assign a plotID type to each row
all_data['geom_str'] = all_data.geometry.apply(lambda x: str(x))
# Step 2: Generate unique IDs for each unique geometry string
all_data['plotId'], unique = pd.factorize(all_data['geom_str'])
# Optional: Drop the temporary 'geom_str' column if you no longer need it
all_data.drop(columns=['geom_str'], inplace=True)
# define df as a pd = to all_data keeping only columns contaiintg either DSM. topoDiv or R_B, plot_type or year
return all_data
# Specify the directory to search for .gpkg files, replace 'your_directory_path' with your actual directory path
version = 8
project_path = r"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline"
# v2: inputs and outputs inside the update folder (the parent of PB Codes); Full_dataset.gpkg is
# Relocated_plots_2\benchmark_final_v9.gpkg (same sha1)
update_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
directory_path = os.path.join(update_dir, 'unchanged', 'Full_dataset.gpkg')
out_path = os.path.join(update_dir, 'run_outputs')
mapbiomas_path = os.path.join(update_dir, 'mapbiomas2018_clip.tif')
all_data = gpd.read_file(directory_path)
# remove last band because that is the date
#all_data = all_data.drop(columns=['B6'])
# If there are rows with the same coordinates and the aim is to have a unique identifier for each unique set of coordinates,
# we first need to identify unique points based on their geometry.
# We will create a new column 'plotID' that assigns a unique ID based on the unique combination of coordinates.
# First, convert the geometry to a string representation to facilitate identifying unique geometries
all_data['geom_str'] = all_data['geometry'].apply(lambda x: str(x))
# Now, assign a unique ID to each unique geometry
unique_geom = {geom: idx + 1 for idx, geom in enumerate(all_data['geom_str'].unique())}
all_data['plotId'] = all_data['geom_str'].map(unique_geom)
# Dropping the temporary 'geom_str' column
all_data.drop(columns=['geom_str'], inplace=True)
#remove geometry
df = all_data.drop(columns=['geometry'])
# create a table with the number of Nan per year divided by plot_type
df_nan = df.isna().groupby([df.Year, df.plot_type]).sum()
df_nan.to_csv(fr'{out_path}/df_nan_v{version}.csv')
# manually select the years for which we have the best coverage of data, based off df_nan
years = [2008, 2010, 2011, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023]
# keep only data from these years
df = df[df.Year.isin(years)]
# which plotId has at least 1 nan in B0, B1, B2, B3, B4, B5
plotId_nan = df[df[['B0', 'B1', 'B2', 'B3', 'B4', 'B5']].isna().any(axis=1)].plotId.unique()
# remove these plotId
df = df[~df.plotId.isin(plotId_nan)]
# remove B0 and B6 columns
df = df.drop(columns=[ 'B6'])
# count nan in each column
print(df.isna().sum())
# remove rows that have NaN values anywhere but in AGB
df = df.dropna(subset=[col for col in df.columns if col != 'AGB'])
df['Year'] = df['Year'].astype(int)
# selecet 45 project and 450 control
# df = select_project_plots(df, n_plots=45)
# get the subset of datapoints with AGB, and calculate the endmembers
df_agb = df.copy().dropna()
df_agb = df_agb[['plotId', 'Year', 'AGB', 'B0', 'B1', 'B2', 'B3', 'B4', 'B5']]
# remove AGB < 1
df_agb = df_agb[df_agb.AGB > np.quantile(df_agb.AGB, .025)]
df_agb = df_agb[df_agb.AGB < np.quantile(df_agb.AGB, .975)]
# v2: Souza NDFI with the fixed CODED endmembers, in place of the NMF endmembers and the band combination chosen by
# correlation with the AGB
df_agb['SI'] = souza_ndfi(df_agb[['B0', 'B1', 'B2', 'B3', 'B4', 'B5']].values)
df_agb = df_agb[np.isfinite(df_agb['SI'])]
# save to file
# keep only plotID, Year, AGB and SI
df_agb_ = df_agb[['plotId', 'Year', 'AGB', 'SI']]
df_agb_.to_csv(fr'{out_path}/df_agb_v{version}.csv', index=False)
# pick the landcover of the values with within plot type
landcover = df[df.plot_type == 'within'].landcover_value.unique()
df = df[df.landcover_value.isin(landcover)]
# v2: MapBiomas Colombia Collection 3.0, 2018: donors keep only the classes held by at least 5 % of the project
# plots, on top of the screen above
with rasterio.open(mapbiomas_path) as src:
mapbiomas_band = src.read(1)
plot_points = all_data.drop_duplicates('plotId').to_crs(src.crs)
rows, cols = rasterio.transform.rowcol(src.transform, plot_points.geometry.x.values, plot_points.geometry.y.values)
mapbiomas = pd.Series(mapbiomas_band[np.asarray(rows), np.asarray(cols)], index=plot_points['plotId'].values)
project_classes = df[df.plot_type == 'within'].drop_duplicates('plotId')['plotId'].map(mapbiomas).value_counts()
dominant_classes = project_classes[project_classes >= 0.05 * project_classes.sum()].index
print(f'v2: MapBiomas 2018 classes of the project plots {project_classes.to_dict()}; dominant {list(dominant_classes)}')
df = df[(df.plot_type == 'within') | df['plotId'].map(mapbiomas).isin(dominant_classes)]
df = df.drop(columns=['AGB'])
df = df.dropna().reset_index(drop=True)
# get unique landuse types
i = df.landcover_value.unique()
# keep only rows with landcover_value in i
df_landcover = df.copy()[df.landcover_value.isin(i)]
# keep only rows with Year in [2014, 2015, 2016, 2017, 2018, 2019]
years = [2008, 2011, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023]
#years = [2008, 2014, 2017, 2018, 2019, 2020, 2021, 2022, 2023]
df_landcover = df_landcover[df_landcover.Year.isin(years)]
# reset index of the dataframe
df_landcover = df_landcover.reset_index(drop=True)
# fill plotId with the index
# count the number of unique plotId in within
print(df_landcover[df_landcover.plot_type == 'within'].plotId.nunique())
# count unique plotId with plot_type == within
print(df_landcover[df_landcover.plot_type == 'outside'].plotId.nunique())
#save df_landcover to gpkg
tmp_gdf = all_data[all_data.Point_ID.isin(df_landcover.Point_ID)]
# keep only the years
tmp_gdf = tmp_gdf[tmp_gdf.Year.isin(years)]
# count unique plotId with plot_type == within in tmp_gdf
print(tmp_gdf[tmp_gdf.plot_type == 'within'].plotId.nunique())
tmp_gdf.to_file(fr'{out_path}/df_no_nan_v{version}.gpkg', driver='GPKG')
# select 30 plots randomly from the plot_type ==within
historic_threshold = 2018
gpkg_dir = out_path # v2
# v2: no NMF endmembers or band combination (fixed CODED endmembers inside analyze)
final_table, final_table_avg, delta_si_control_history, p_value_control_history, std_err_control_history, delta_si_control_mean_history,p_value_control_mean_history, std_err_control_mean_history = analyze(df_landcover, historic_threshold=historic_threshold, endmember_spectra =None,
best_combination=None, n_controls_per_plot = 4,
out_path=out_path, filename=fr'{out_path}/i_summary_v{version}', get_the_smallest_distance=45, gpkg_dir=gpkg_dir, source_gpkg=tmp_gdf)
print(final_table)
print(delta_si_control_history)
print(p_value_control_history)