Now we need to take the sparse fishnets and extract our remote sensed data. to do so we will do the following:

import all the

In [1]:

import ee
import pandas as pd
import geopandas as gpd
import geemap
import os
from tqdm import tqdm
from shapely import wkt
from shapely.geometry import mapping

#ee.Authenticate()
ee.Initialize()

In [4]:

os.chdir(os.path.abspath(os.path.join(os.getcwd(), os.pardir)))

In [5]:

os.getcwd()

'/home/datascience/herbivore_nutrient_interactions'

In [6]:

# ISRIC SoilGrid

# Initialize the Earth Engine library
ee.Initialize()

# Define the list of asset paths and corresponding soil property names
soil_properties = {
    'bdod_mean': 'bdod',
    'cec_mean': 'cec',
    'cfvo_mean': 'cfvo',
    'clay_mean': 'clay',
    'nitrogen_mean': 'nitrogen',
    'ocd_mean': 'ocd',
    'ocs_mean': 'ocs',
    'phh2o_mean': 'phh2o',
    'sand_mean': 'sand',
    'soc_mean': 'soc',
    'silt_mean': 'silt'
}

# Define the depths of interest
depths = ['0-5cm', '5-15cm']

def median_first_two_depths(asset_path, property_name):
    # Load the image
    image = ee.Image(asset_path)
    
    # Handle the special case for 'ocs_mean' which has only one band
    if property_name == 'ocs':
        # Return the image directly
        median_image = image.rename([f'{property_name}_single_band'])
    else:
        # Construct the band names for other properties
        band_names = [f'{property_name}_{depth}_mean' for depth in depths]
        
        # Select the bands of interest
        selected_bands = image.select(band_names)
        
        # Compute the median of the selected bands
        median_image = selected_bands.reduce(ee.Reducer.median())
        
        # Rename the band to include the soil property name
        median_image = median_image.rename([f'{property_name}_median'])
    
    return median_image

# Create a list to hold the median images
median_images = [median_first_two_depths(f'projects/soilgrids-isric/{key}', value) for key, value in soil_properties.items()]

# Combine the median images into an ImageCollection
isric_soil_grid = ee.ImageCollection(median_images).toBands()

In [7]:

# Australia soil grid


# Load the Australian Soil Grid image collection.
soil_collection = ee.ImageCollection("CSIRO/SLGA")

# Define the depth ranges and prefixes
depth_ranges = ['000_005', '005_015']
prefixes = ['SOC', 'AWC', 'BDW', 'CLY', 'ECE', 'NTO', 'pHc', 'PTO', 'SLT', 'SND']

# Initialize an empty list to store the band names
selected_images_redundant = []
selected_images_final = []
# Loop through each prefix and depth range to construct the band names
for prefix in prefixes:
    for depth_range in depth_ranges:
        selected_images_final.append(f'{prefix}_{depth_range}_EV')
        selected_images_redundant.append(f'{prefix}_{prefix}_{depth_range}_EV')


# take the imageCollection
# convert to iamge with bands
# select and rename the bands

slga_image = soil_collection.toBands().select(selected_images_redundant)

In [8]:

# Bioclim variables
bioclim = ee.Image("WORLDCLIM/V1/BIO")

# Ecoregion
ecoregion = ee.FeatureCollection("RESOLVE/ECOREGIONS/2017").select('BIOME_NAME','ECO_NAME','ECO_ID','BIOME_NUM')

# Elevation
elevation = ee.Image("USGS/SRTMGL1_003")

# woody vegetation pixel cover

woody_veg = ee.ImageCollection("NASA/MEASURES/GFCC/TC/v3").select('tree_canopy_cover').median()

In [16]:

ecoregions = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017')

In [9]:

first = isric_soil_grid.addBands(elevation)
second = first.addBands(bioclim)
third = second.addBands(slga_image)
final = third.addBands(woody_veg)
final

Now that we have all the bands and features together, lets extract to the grids

lets import both the APL and DL csvs and convert them into geopandas

In [7]:

# create APL dataframe

apl_df = pd.read_csv('data/processed/spatial_modeling/aplc_data_aggregated_to_polygon_grid.csv')
apl_df['geometry'] = apl_df['geometry'].apply(wkt.loads)
apl_gdf = gpd.GeoDataFrame(apl_df, geometry='geometry')

In [11]:


# Function to get median band values for a single polygon
def get_median_values(polygon):
    ee_polygon = ee.Geometry.Polygon(mapping(polygon)['coordinates'])
    median_values = final.reduceRegion(
        reducer=ee.Reducer.median(),
        geometry=ee_polygon,
        scale=30,  # Change this to the appropriate scale for your data
        maxPixels=1e9
    ).getInfo()
    return median_values



# Process the DataFrame in chunks
chunk_size = 50
apl_chunks = [apl_gdf[i:i + chunk_size] for i in range(0, apl_gdf.shape[0], chunk_size)]

apl_results = []

for chunk in tqdm(apl_chunks, desc="Processing chunks"):
    for index, row in chunk.iterrows():
        polygon_id = row['polygon_id']  # Assuming 'polygon_id' is the column name
        polygon = row['geometry']
        median_values = get_median_values(polygon)
        median_values['geometry'] = polygon
        median_values['polygon_id'] = polygon_id  # Add polygon_id to the results
        apl_results.append(median_values)

Processing chunks: 100%|██████████| 1984/1984 [3:11:44<00:00,  5.80s/it]

In [19]:

import ee
from shapely.geometry import mapping
from tqdm import tqdm

# Initialize Earth Engine API
ee.Initialize()

# Load the ecoregions FeatureCollection
ecoregions = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017')

# Function to get median band values and ecoregion name for a single polygon
def get_median_values(polygon):
    ee_polygon = ee.Geometry.Polygon(mapping(polygon)['coordinates'])
    
    # Get median band values from the final dataset
    median_values = final.reduceRegion(
        reducer=ee.Reducer.median(),
        geometry=ee_polygon,
        scale=30,  # Adjust this scale as necessary
        maxPixels=1e9
    ).getInfo()
    
    # Get the first ecoregion that intersects with the polygon
    intersecting_ecoregion = ecoregions.filterBounds(ee_polygon).first()
    
    # Retrieve the 'ECO_NAME' property (or any relevant property) from the intersecting ecoregion
    if intersecting_ecoregion:
        ecoregion_name = intersecting_ecoregion.get('ECO_NAME').getInfo()
    else:
        ecoregion_name = 'Unknown'  # If no intersection found, use 'Unknown'
    
    # Add the ecoregion name to the median values
    median_values['ecoregion'] = ecoregion_name
    
    return median_values

# Process the DataFrame in chunks
chunk_size = 100
apl_chunks = [apl_gdf[i:i + chunk_size] for i in range(0, apl_gdf.shape[0], chunk_size)]

apl_results = []

for chunk in tqdm(apl_chunks, desc="Processing chunks"):
    for index, row in chunk.iterrows():
        polygon_id = row['polygon_id']  # Assuming 'polygon_id' is the column name
        polygon = row['geometry']
        
        # Get median values and ecoregion name for the current polygon
        median_values = get_median_values(polygon)
        
        # Add additional information to the results
        median_values['geometry'] = polygon
        median_values['polygon_id'] = polygon_id
        
        apl_results.append(median_values)

# apl_results now contains median band values and the corresponding ecoregion name

Processing chunks: 100%|██████████| 992/992 [11:14:20<00:00, 40.79s/it]

In [20]:

apl_results_df = pd.DataFrame(apl_results)
apl_results_df.to_csv('data/processed/spatial_modeling/apl_remote_sense_extraction.csv', index=False)

Now lets relate the two dataframes together and serialize

In [21]:

apl_results_df = pd.read_csv('data/processed/spatial_modeling/apl_remote_sense_extraction.csv')

apl_survey_remote_sense_data = apl_df.merge(apl_results_df, on='polygon_id', how='left')

apl_survey_remote_sense_data.to_csv('data/processed/spatial_modeling/apl_survey_remote_sense_data.csv', index=False)

Now we have the australian plague locust dataset as a dataframe

3-2 spatial modeling - google earth extraction

Extracting GEE data to polygons

Now that we have all the bands and features together, lets extract to the grids

Now lets relate the two dataframes together and serialize