Timeseries extraction

Uncomment the following line to install geeagri if needed.

In [ ]:

# !pip install geeagri

# !pip install geeagri

Import libraries¶

In [ ]:

import ee
import geemap
from geeagri.extract import (
    extract_timeseries_to_point,
    extract_timeseries_to_polygon,
    TimeseriesExtractor,
)

import matplotlib.pyplot as plt
import seaborn as sns

plt.rcParams["font.family"] = "DeJavu Serif"
plt.rcParams["font.serif"] = "Times New Roman"

import ee import geemap from geeagri.extract import ( extract_timeseries_to_point, extract_timeseries_to_polygon, TimeseriesExtractor, ) import matplotlib.pyplot as plt import seaborn as sns plt.rcParams["font.family"] = "DeJavu Serif" plt.rcParams["font.serif"] = "Times New Roman"

Initialize a Map object¶

Authenticate and initialize Earth Engine. If it doesn't work, specify a project name

In [ ]:

# ee.Authenticate()
# ee.Initialize(project='your-project-id')

Map = geemap.Map()
Map

# ee.Authenticate() # ee.Initialize(project='your-project-id') Map = geemap.Map() Map

Import region of interest¶

In [ ]:

bbox = [-100.612793, 29.084977, -95.679932, 31.896214]
polygon = ee.Geometry.BBox(*bbox)
polygon_style = {"color": "red", "width": 1}
Map.addLayer(polygon, polygon_style, "Polygon")
Map.centerObject(polygon, 8)

lon, lat = -98.15, 30.50
point = ee.Geometry.Point([lon, lat])
Map.addLayer(point, {"color": "blue"}, "Point")

bbox = [-100.612793, 29.084977, -95.679932, 31.896214] polygon = ee.Geometry.BBox(*bbox) polygon_style = {"color": "red", "width": 1} Map.addLayer(polygon, polygon_style, "Polygon") Map.centerObject(polygon, 8) lon, lat = -98.15, 30.50 point = ee.Geometry.Point([lon, lat]) Map.addLayer(point, {"color": "blue"}, "Point")

Parallel export of timeseries for large sample sets with TimeseriesExtractor¶

In [ ]:

era5_land = ee.ImageCollection("ECMWF/ERA5_LAND/DAILY_AGGR")

# Generate or load sample points
sample = (
    era5_land.first()
    .select("temperature_2m_min")
    .sample(
        region=polygon,
        scale=11132,  # ~11 km
        numPixels=100,  # number of samples
        seed=42,
        dropNulls=True,
        tileScale=16,
        geometries=True,
    )
)

# Convert samples to GeoDataFrame
sample_gdf = geemap.ee_to_gdf(sample)
sample_gdf["ID"] = sample_gdf.index
sample_gdf = sample_gdf[["ID", "geometry"]]

print(f"Sample size: {sample_gdf.shape[0]}")
sample_gdf.head()

era5_land = ee.ImageCollection("ECMWF/ERA5_LAND/DAILY_AGGR") # Generate or load sample points sample = ( era5_land.first() .select("temperature_2m_min") .sample( region=polygon, scale=11132, # ~11 km numPixels=100, # number of samples seed=42, dropNulls=True, tileScale=16, geometries=True, ) ) # Convert samples to GeoDataFrame sample_gdf = geemap.ee_to_gdf(sample) sample_gdf["ID"] = sample_gdf.index sample_gdf = sample_gdf[["ID", "geometry"]] print(f"Sample size: {sample_gdf.shape[0]}") sample_gdf.head()

In [ ]:

# Extract timeseries in parallel for all samples
ts_extractor = TimeseriesExtractor(
    image_collection=era5_land,
    sample_gdf=sample_gdf,
    identifier="ID",
    out_dir="test",  # output directory
    selectors=[
        "temperature_2m_min",
        "temperature_2m_max",
        "total_precipitation_sum",
        "surface_solar_radiation_downwards_sum",
    ],
    scale=11132,
    num_processes=20,  # parallel processes
    start_date="2000-01-01",
    end_date="2010-01-01",
)

# Run extraction
ts_extractor.extract_timeseries()

# Extract timeseries in parallel for all samples ts_extractor = TimeseriesExtractor( image_collection=era5_land, sample_gdf=sample_gdf, identifier="ID", out_dir="test", # output directory selectors=[ "temperature_2m_min", "temperature_2m_max", "total_precipitation_sum", "surface_solar_radiation_downwards_sum", ], scale=11132, num_processes=20, # parallel processes start_date="2000-01-01", end_date="2010-01-01", ) # Run extraction ts_extractor.extract_timeseries()

Export single timeseries from climate data (ERA5-Land Daily)¶

In [ ]:

era5_land = ee.ImageCollection("ECMWF/ERA5_LAND/DAILY_AGGR")

era5_land_point_ts = extract_timeseries_to_point(
    lat=lat,
    lon=lon,
    image_collection=era5_land,
    start_date="2010-01-01",
    end_date="2015-01-01",
    band_names=[
        "temperature_2m_min",
        "temperature_2m_max",
        "total_precipitation_sum",
        "surface_solar_radiation_downwards_sum",
    ],
    scale=11132,
)

era5_land_polygon_ts = extract_timeseries_to_polygon(
    polygon=polygon,
    image_collection=era5_land,
    start_date="2010-01-01",
    end_date="2015-01-01",
    band_names=[
        "temperature_2m_min",
        "temperature_2m_max",
        "total_precipitation_sum",
        "surface_solar_radiation_downwards_sum",
    ],
    scale=11132,
    reducer="MEAN",
)

era5_land = ee.ImageCollection("ECMWF/ERA5_LAND/DAILY_AGGR") era5_land_point_ts = extract_timeseries_to_point( lat=lat, lon=lon, image_collection=era5_land, start_date="2010-01-01", end_date="2015-01-01", band_names=[ "temperature_2m_min", "temperature_2m_max", "total_precipitation_sum", "surface_solar_radiation_downwards_sum", ], scale=11132, ) era5_land_polygon_ts = extract_timeseries_to_polygon( polygon=polygon, image_collection=era5_land, start_date="2010-01-01", end_date="2015-01-01", band_names=[ "temperature_2m_min", "temperature_2m_max", "total_precipitation_sum", "surface_solar_radiation_downwards_sum", ], scale=11132, reducer="MEAN", )

In [ ]:

# Plot the data
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 6))
axes = axes.flatten()

sns.lineplot(
    data=era5_land_point_ts,
    x="time",
    y="temperature_2m_max",
    c="r",
    linewidth=0.5,
    ax=axes[0],
    label="Tmax",
)
sns.lineplot(
    data=era5_land_point_ts,
    x="time",
    y="temperature_2m_min",
    c="b",
    linewidth=0.5,
    ax=axes[0],
    label="Tmax",
)
axes[0].set_ylabel("Values")
axes[0].legend()
axes[0].set_title("Daily timeseries of climate data based on point")

sns.lineplot(
    data=era5_land_polygon_ts,
    x="time",
    y="temperature_2m_max",
    c="r",
    linewidth=0.5,
    ax=axes[1],
    label="Tmax",
)
sns.lineplot(
    data=era5_land_polygon_ts,
    x="time",
    y="temperature_2m_min",
    c="b",
    linewidth=0.5,
    ax=axes[1],
    label="Tmax",
)
axes[1].set_ylabel("Values")
axes[1].legend()
axes[1].set_title("Daily timeseries of climate data based on polygon")

plt.tight_layout()
plt.show()

# Plot the data fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 6)) axes = axes.flatten() sns.lineplot( data=era5_land_point_ts, x="time", y="temperature_2m_max", c="r", linewidth=0.5, ax=axes[0], label="Tmax", ) sns.lineplot( data=era5_land_point_ts, x="time", y="temperature_2m_min", c="b", linewidth=0.5, ax=axes[0], label="Tmax", ) axes[0].set_ylabel("Values") axes[0].legend() axes[0].set_title("Daily timeseries of climate data based on point") sns.lineplot( data=era5_land_polygon_ts, x="time", y="temperature_2m_max", c="r", linewidth=0.5, ax=axes[1], label="Tmax", ) sns.lineplot( data=era5_land_polygon_ts, x="time", y="temperature_2m_min", c="b", linewidth=0.5, ax=axes[1], label="Tmax", ) axes[1].set_ylabel("Values") axes[1].legend() axes[1].set_title("Daily timeseries of climate data based on polygon") plt.tight_layout() plt.show()

Export single timeseries from NDVI data (MOD13Q1.061 Terra Vegetation Indices 16-Day)¶

In [ ]:

modis_ndvi = ee.ImageCollection("MODIS/061/MOD13Q1")

modis_ndvi_point_ts = extract_timeseries_to_point(
    lat=lat,
    lon=lon,
    image_collection=modis_ndvi,
    start_date="2010-01-01",
    end_date="2015-01-01",
    band_names=["NDVI", "EVI"],
    scale=250,
)

modis_ndvi_polygon_ts = extract_timeseries_to_polygon(
    polygon=polygon,
    image_collection=modis_ndvi,
    start_date="2010-01-01",
    end_date="2015-01-01",
    band_names=["NDVI", "EVI"],
    scale=250,
    reducer="MEAN",
)

# Apply scale factors
modis_ndvi_point_ts[["NDVI", "EVI"]] = modis_ndvi_point_ts[["NDVI", "EVI"]] * 0.0001
modis_ndvi_polygon_ts[["NDVI", "EVI"]] = modis_ndvi_polygon_ts[["NDVI", "EVI"]] * 0.0001

modis_ndvi = ee.ImageCollection("MODIS/061/MOD13Q1") modis_ndvi_point_ts = extract_timeseries_to_point( lat=lat, lon=lon, image_collection=modis_ndvi, start_date="2010-01-01", end_date="2015-01-01", band_names=["NDVI", "EVI"], scale=250, ) modis_ndvi_polygon_ts = extract_timeseries_to_polygon( polygon=polygon, image_collection=modis_ndvi, start_date="2010-01-01", end_date="2015-01-01", band_names=["NDVI", "EVI"], scale=250, reducer="MEAN", ) # Apply scale factors modis_ndvi_point_ts[["NDVI", "EVI"]] = modis_ndvi_point_ts[["NDVI", "EVI"]] * 0.0001 modis_ndvi_polygon_ts[["NDVI", "EVI"]] = modis_ndvi_polygon_ts[["NDVI", "EVI"]] * 0.0001

In [ ]:

# Plot the data
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 6))
axes = axes.flatten()

sns.lineplot(
    data=modis_ndvi_point_ts,
    x="time",
    y="NDVI",
    c="green",
    marker="o",
    markersize=5,
    linewidth=1,
    ax=axes[0],
    label="NDVI",
)
sns.lineplot(
    data=modis_ndvi_point_ts,
    x="time",
    y="EVI",
    c="orange",
    marker="o",
    markersize=5,
    linewidth=1,
    ax=axes[0],
    label="EVI",
)
axes[0].set_ylabel("Values")
axes[0].legend()
axes[0].set_title("Daily timeseries of NDVI and EVI data based on point")

sns.lineplot(
    data=modis_ndvi_polygon_ts,
    x="time",
    y="NDVI",
    c="green",
    marker="o",
    markersize=5,
    linewidth=1,
    ax=axes[1],
    label="NDVI",
)
sns.lineplot(
    data=modis_ndvi_polygon_ts,
    x="time",
    y="EVI",
    c="orange",
    marker="o",
    markersize=5,
    linewidth=1,
    ax=axes[1],
    label="EVI",
)
axes[1].set_ylabel("Values")
axes[1].legend()
axes[1].set_title("Daily timeseries of NDVI and EVI data based on polygon")

plt.tight_layout()
plt.show()

# Plot the data fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(12, 6)) axes = axes.flatten() sns.lineplot( data=modis_ndvi_point_ts, x="time", y="NDVI", c="green", marker="o", markersize=5, linewidth=1, ax=axes[0], label="NDVI", ) sns.lineplot( data=modis_ndvi_point_ts, x="time", y="EVI", c="orange", marker="o", markersize=5, linewidth=1, ax=axes[0], label="EVI", ) axes[0].set_ylabel("Values") axes[0].legend() axes[0].set_title("Daily timeseries of NDVI and EVI data based on point") sns.lineplot( data=modis_ndvi_polygon_ts, x="time", y="NDVI", c="green", marker="o", markersize=5, linewidth=1, ax=axes[1], label="NDVI", ) sns.lineplot( data=modis_ndvi_polygon_ts, x="time", y="EVI", c="orange", marker="o", markersize=5, linewidth=1, ax=axes[1], label="EVI", ) axes[1].set_ylabel("Values") axes[1].legend() axes[1].set_title("Daily timeseries of NDVI and EVI data based on polygon") plt.tight_layout() plt.show()