Script 1 — Data retriever
Generates the plots and assembles the dataset on which everything else runs: the random draw within the donor pool area, the Landsat reflectance by band and year, the biomass and the land cover.
The listing below is vm0047_data_retriever.py exactly as delivered: 388 lines, unmodified. It is read from the file at render time, so the document and the script cannot drift apart.
⬇ Download vm0047_data_retriever.py
#(c) Terra Global Capital CONFIDENTIAL
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
import random
from geopandas import GeoDataFrame
from shapely.geometry import Point, box
import geopandas as gpd
import rasterio
from rasterio.features import geometry_mask
from rasterio.transform import from_origin
import rasterio
import geopandas as gpd
import rasterio
import numpy as np
from rasterio.transform import from_origin
from rasterio.enums import Resampling
import geopandas as gpd
import numpy as np
import os
#import ee
def generate_random_points_within_extent(extent, num_points, shapefile):
minx, miny, maxx, maxy = extent
points = []
while len(points) < num_points:
x = random.uniform(minx, maxx)
y = random.uniform(miny, maxy)
point = Point(x, y)
if shapefile.geometry.contains(point).any():
points.append(point)
return points
def generate_random_points_around_centroid(centroid, extent, num_points, distance_km):
buffered_area = centroid.buffer(distance_km * 1000) # rough conversion to degrees
extent_box = box(*extent)
outside_area = buffered_area.difference(extent_box)
points = []
while len(points) < num_points:
x = random.uniform(outside_area.bounds[0], outside_area.bounds[2])
y = random.uniform(outside_area.bounds[1], outside_area.bounds[3])
point = Point(x, y)
if outside_area.contains(point):
points.append(point)
return points
def get_plots_and_controls(shapefile_path, landcover_path, output_path, radius = 50, plots_pool = 30, epsg = 32619, control_per_plot = 100):
# Get centroid and extent
shapefile = gpd.read_file(shapefile_path)
centroid = shapefile.geometry.centroid.to_crs(epsg=epsg).iloc[0] # Assuming WGS84 projection
extent = shapefile.geometry.total_bounds
# Generate 30 random locations within the shapefile
points_within = generate_random_points_within_extent(extent, plots_pool, shapefile)
# Generate 150 random locations within 100 km but outside the shapefile
points_outside = generate_random_points_around_centroid(centroid, extent, plots_pool*control_per_plot, radius)
# Convert the list of Points to a GeoSeries
points_within_geo = gpd.GeoSeries(points_within)
points_outside_geo = gpd.GeoSeries(points_outside)
# Create GeoDataFrames
gdf_within = GeoDataFrame(geometry=points_within_geo)
gdf_outside = GeoDataFrame(geometry=points_outside_geo)
# Add "plot_type" column
gdf_within['plot_type'] = 'within'
gdf_outside['plot_type'] = 'outside'
# Optionally, combine the two DataFrames
combined_gdf = pd.concat([gdf_within, gdf_outside]).reset_index(drop=True)
# add crs to the combined_gdf
combined_gdf.crs = shapefile.crs
# Load the raster from 'landcover_path'
with rasterio.open(landcover_path) as src:
raster_crs = src.crs
raster_transform = src.transform
# Reproject the GeoDataFrame to match the raster CRS
if combined_gdf.crs != raster_crs:
combined_gdf = combined_gdf.to_crs(raster_crs)
# Initialize an empty list to store raster values
raster_values = []
# Extract raster values at each point's coordinates
for idx, row in combined_gdf.iterrows():
# Convert the point geometry to raster space
row_x, row_y = row.geometry.x, row.geometry.y
col, row = ~raster_transform * (row_x, row_y)
col, row = int(col), int(row)
# Read the raster value at the point location
for val in src.sample([(row_x, row_y)]):
raster_values.append(val[0] if src.count == 1 else val) # Handle single or multi-band rasters
# Append the raster values to the GeoDataFrame
combined_gdf['landcover_value'] = raster_values
#
# Example: Save to a new shapefile
combined_gdf.to_file(output_path, driver="GPKG")
return combined_gdf
def is_collection(data_product):
"""
Check if the given Earth Engine data product is an ImageCollection.
Parameters:
data_product (str): Earth Engine data product ID.
Returns:
bool: True if the data product is an ImageCollection, False otherwise.
"""
try:
# Attempt to use a method unique to ImageCollection
ee.ImageCollection(data_product).filterDate('2000-01-01', '2000-01-02')
return True
except:
return False
def flatten_band_info(band_info_list):
"""Flatten a list of dictionaries into a single dictionary."""
flattened_info = {}
for item in band_info_list:
flattened_info.update(item)
return flattened_info
def maskL(image):
qa = image.select('QA_PIXEL')
# Use appropriate bits for cloud and shadow masking
cloud_shadow = qa.bitwiseAnd(1 << 3)
clouds = qa.bitwiseAnd(1 << 5)
mask = cloud_shadow.Or(clouds).eq(0)
return image.updateMask(mask)
def maskClouds(image):
fmask = image.select('Fmask')
mask = fmask.bitwiseAnd(1 << 1).eq(0) \
.And(fmask.bitwiseAnd(1 << 2).eq(0)) \
.And(fmask.bitwiseAnd(1 << 3).eq(0))
return image.updateMask(mask)
def process_data_product(data_product, feature, start_date, end_date):
# Initialize an empty dictionary to store band information
band_info = {}
# Define AOI and point
point = ee.Geometry.Point([feature.geometry.x, feature.geometry.y])
# define the AOI as point buffered 10 pixels
AOI = point.buffer(10)
# Check if data_product is a collection or an image
if 'LANDSAT' in data_product:
collection = ee.ImageCollection(data_product).filterDate(start_date, end_date).map(maskL)
image = collection.median()
band_info = image.sample(point, 30).first()
if band_info.getInfo(): # This checks if the sampled feature is not None
band_info = band_info.toDictionary().getInfo()
# Add all bands starting with 'SR_' to the band_info dictionary
for key, value in band_info.items():
if 'SR_' in key:
band_info[key] = value
else:
# If no data is found, set band_info to None
band_info = None # Handle the case where no data is found
elif 'HLS' in data_product:
collection = ee.ImageCollection('NASA/HLS/HLSL30/v002').filterDate(start_date, end_date).filterBounds(AOI).map(maskClouds)
image = collection.median()
band_info = image.sample(point, 30).first()
if band_info.getInfo(): # This checks if the sampled feature is not None
band_info = band_info.toDictionary().getInfo()
# Add all bands starting with 'SR_' to the band_info dictionary
for key, value in band_info.items():
if 'B' in key:
band_info[key] = value
elif 'AW3D30' in data_product:
dataset = ee.Image(data_product)
band_info = dataset.sample(point, 30).first()
if band_info.getInfo(): # This checks if the sampled feature is not None
band_info = band_info.toDictionary().getInfo()
else:
# If no data is found, set band_info to None
band_info = None # Handle the case where no data is found
#band_info = flatten_band_info(band_info)
elif 'topoDiversity' in data_product:
dataset = ee.Image(data_product)
band_info = dataset.sample(point, 30).first()
if band_info.getInfo(): # This checks if the sampled feature is not None
band_info = band_info.toDictionary().getInfo()
else:
# If no data is found, set band_info to None
band_info = None # Handle the case where no data is found
elif 'Landcover' in data_product:
dataset = ee.Image(data_product)
band_info = dataset.sample(point, 30).first()
if band_info.getInfo(): # This checks if the sampled feature is not None
band_info = band_info.toDictionary().getInfo()
else:
# If no data is found, set band_info to None
band_info = None # Handle the case where no data is found
return band_info
def build_dataset_for_benchmark(shapefile_path, landcover_path, output_path, project_path, version = 0, plots_pool = 50, control_per_plot = 100, epsg = 32619, project_id = 'earth-396713', first_year = 2008, last_year = 2023, reflectance_from_sr = True, source_path = r"\1_Source", agb_path = fr"\1_Source\agb", refine_boundaries = None, reuse_data = None):
if reuse_data is not None:
final_gpd = gpd.read_file(reuse_data)
else:
final_gpd = get_plots_and_controls(shapefile_path, landcover_path, output_path, epsg = epsg, plots_pool = plots_pool, control_per_plot = control_per_plot)
#final_gpd = gpd.read_file(r"C:\Users\sergio.marconi\Documents\GitHub\dynamic_baseline\benchmark_v2\benchmark_2010_v0.gpkg")
# keep only geometry plot_type landcover_value columns
final_gpd = final_gpd[['geometry', 'plot_type', 'landcover_value']]
years = range(first_year, last_year)
# Load Geopandas DataFrame
final_gpd = final_gpd.to_crs(epsg=4326) # Reproject to WGS 84
final_gpd["plotId"] = final_gpd.index
# List of GEE data products
if refine_boundaries is not None:
for conditions_shp in refine_boundaries:
conditions = gpd.read_file(conditions_shp)
conditions = conditions.to_crs(epsg=4326)
# remove from final_gpd all points that are not within the conditions region
tmp_ = gpd.sjoin(final_gpd, conditions, how="inner", predicate='within')
# filter out plotId that don't show in tmp_.plotId
final_gpd = final_gpd[final_gpd.plotId.isin(tmp_.plotId)]
if reflectance_from_sr:
# get the files from 1_Source
source_path = source_path
# list files in the folder
files = os.listdir(source_path)
# order the files by date
files.sort()
# Initialize an empty DataFrame for the long format dataset
long_format_df = pd.DataFrame(columns=['Point_ID', 'Year', 'Band', 'Value'])
# Initialize a list to collect data
data_list = []
for fl in files:
if '.tif' in fl:
# Get the year from the file name
year = int(fl.split('_')[1][:-10]) # Adjust slicing if needed
with rasterio.open(source_path + '/' + fl) as src:
for idx, row in final_gpd.iterrows():
# Convert the point geometry to raster space
row_x, row_y = row.geometry.x, row.geometry.y
col, row = ~src.transform * (row_x, row_y)
col, row = int(col), int(row)
if (0 <= col < src.width) and (0 <= row < src.height):
values = src.read(window=rasterio.windows.Window(col, row, 1, 1), masked=True).flatten()
for band_idx, value in enumerate(values):
if value is np.ma.masked:
# Handle missing data if necessary
value = np.nan
# Collect data in a list
data_list.append({
'Point_ID': idx, # Adjust if there's a better identifier
'Year': year,
'Band': f"B{band_idx}",
'Value': value
})
else:
# Optionally handle points outside raster bounds
pass
# Create a DataFrame from the collected data
long_format_df = pd.DataFrame(data_list)
# list the unique Point_ID that show at least 1 nan
points_with_nan = long_format_df[long_format_df.Value.isna()].Point_ID.unique()
# remove points with nan
#long_format_df = long_format_df[~long_format_df.Point_ID.isin(points_with_nan)]
# Ensure correct data types, if necessary
long_format_df['Year'] = long_format_df['Year'].astype(int)
long_format_df['Value'] = long_format_df['Value'].astype(float)
# now add the AGB
source_path = agb_path
# list files in the folder
files = os.listdir(source_path)
# order the files by date
# Initialize an empty DataFrame for the long format dataset
long_format_agb = pd.DataFrame(columns=['Point_ID', 'Year', 'Band', 'Value'])
# Initialize a list to collect data
data_list = []
for fl in files:
if '.tif' in fl:
if '.aux.xml' in fl:
continue
# Get the year from the file name
year = int(fl.split('_')[1][-8:-4]) # Adjust slicing if needed
with rasterio.open(source_path + '/' + fl) as src:
for idx, row in final_gpd.iterrows():
# Convert the point geometry to raster space
row_x, row_y = row.geometry.x, row.geometry.y
col, row = ~src.transform * (row_x, row_y)
col, row = int(col), int(row)
if (0 <= col < src.width) and (0 <= row < src.height):
values = src.read(window=rasterio.windows.Window(col, row, 1, 1), masked=True).flatten()
for band_idx, value in enumerate(values):
if value is np.ma.masked:
# Handle missing data if necessary
value = np.nan
# Collect data in a list
data_list.append({
'Point_ID': idx, # Adjust if there's a better identifier
'Year': year,
'Band': f"AGB",
'Value': value
})
else:
# Optionally handle points outside raster bounds
pass
# Create a DataFrame from the collected data
long_format_agb = pd.DataFrame(data_list)
# Ensure correct data types, if necessary
long_format_agb['Year'] = long_format_agb['Year'].astype(int)
long_format_agb['Value'] = long_format_agb['Value'].astype(float)
#append the long_format_agb to the long_format_df
long_format = pd.concat([long_format_df, long_format_agb])
df_wide = long_format.pivot_table(index=["Point_ID", "Year"], columns="Band", values="Value").reset_index()
# Setting Point_ID as the index for the DataFrame to merge
df_wide.set_index('Point_ID', inplace=True)
df_wide['Point_ID'] = df_wide.index
# Merging the wide DataFrame with the GeoDataFrame based on the index
gdf_merged = final_gpd.merge(df_wide, left_index=True, right_index=True, how='left')
#gdf_merged = gdf_merged.dropna()
target_lulcs = gdf_merged[gdf_merged.plot_type == 'within'].landcover_value.unique()
# get unique Point_ID and landcover_value from the gdf_merged and count
# remove all rows where landcover_value is not in target_lulcs
gdf_merged = gdf_merged[gdf_merged.landcover_value.isin(target_lulcs)]
# get unique Point_ID, plot_type and landcover_value from the gdf_merged
check = gdf_merged[['Point_ID', 'plot_type', 'landcover_value']].drop_duplicates()
# count plot_type and landcover_value combinations
check.groupby(['plot_type', 'landcover_value']).size().reset_index(name='counts')
# check the number of unique Point_ID that are within
check = gdf_merged[['Point_ID', 'plot_type']].drop_duplicates()
check.groupby(['plot_type']).size().reset_index(name='counts')
# save to a new file
gdf_merged.to_file(rf"{project_path}\benchmark_final_v{version}.gpkg", driver="GPKG")
#
# Load the shapefile
project_path = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline"
shapefile_path = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline\conditions\project_area.gpkg" # Replace with your shapefile path
landcover_path = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline\1_Source\landuse_history.tif"
output_path = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline\Relocated_plots_2\combined_points.gpkg"
source_landsat = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\1_Source\1_Multispectral\Landsat_PB"
source_agb = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\1_Source\5_Biophysical\AGB"
refine_boundaries_shps = [fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline\Relocated_plots_2\Region_pool_v4.gpkg"]
first_year = 2008
last_year = 2024
# Initialize Earth Engine
version = 9
# Initialize Earth Engine
project_id = 'ee-tgcwindows11'
crs = 'EPSG:4326' # WGS 84 coordinate system
path_dataset = fr"A:\TGC-1643 12Tree - Colombia - La Paz - ARR w harvest\2_RS\5_Dynamic_Baseline\Relocated_plots_2\datset_v{version}.gpkg"
build_dataset_for_benchmark(shapefile_path, landcover_path, output_path, plots_pool = 150,
epsg = 32619, project_id = project_id, project_path = project_path,
first_year = first_year, last_year = last_year, reflectance_from_sr = True,
version = version, source_path = source_landsat, agb_path = source_agb,
refine_boundaries = refine_boundaries_shps, reuse_data = None)