nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Extract values to point from multiple raster with time property

Pva013

Hi, i have multiple rasters of VCI (Vegetation Condition Index) from 1/12/2019 to 31/12/2022 (from December to June) with 16 days interval reset at 01/01 of each year like this image

I also have a set of points with date like this.

I want to extract values to point but with the exact time period (Ex: VCI raster from 1/1/2020 to 17/1/2020 extract to points have date in that period like 5/1/2020). I tried to create a mosaic dataset then set time property and date field to it but after that i have no idea how to continue. It would be great if anyone have other solution to the problem.

Thanks in advance

Find more posts tagged with

Content Management

Desktop

Feature Services

Accepted answers

AustinAverill

I would say the best way to do this would probably be to iterate through the rasters and read in the date period from the name using python. You could then do a selectByAttribute within the loop and then extract values based on that selection. It would look something like this.

import arcpy
from datetime import datetime

arcpy.env.workspace = "path\\to\\fgdb\\with\\rasters"

points = "path\\to\\points"

rasters = arcpy.ListRasters()

for i in range(len(rasters)):
    while i < len(rasters) - 1:
        start_date_str = raster[i].name.replace("VCI_VCI_", "")
        end_date_str = raster[i+1].name.replace("VCI_VCI_, "")
        start_date_format = "%d_%m_%Y"
        final_date_format = "%Y_%m_%d"
        
        start_date = datetime.strptime(start_date_str, start_date_format).strftime(final_date_format)
        end_date = datetime.strptrime(end_date_str, start_date_format).strftime(final_date_format)

        q_str = f"'ACQ_DATE' >= date '{start_date}' AND 'ACQ_DATE' < date '{end_date}'"
        selection = arcpy.management.SelectLayerByAttribute(points, "NEW_SELECTION", q_str)

        arcpy.sa.ExtractMultiValuesToPoints(selection, rasters[i])

This might not be the perfect example, but should get you going in the right direction.

All comments

MervynLotter

After creating your raster mosaic, the next step would be to create a multidimensional layer using the "Make Multidimensional Raster Layer" GP tool. You will them have a set of multidimensional raster tools available to explore your data and extract the necessary information.

See Multidimensional raster data—ArcGIS Pro | Documentation

AustinAverill

import arcpy
from datetime import datetime

arcpy.env.workspace = "path\\to\\fgdb\\with\\rasters"

points = "path\\to\\points"

rasters = arcpy.ListRasters()

for i in range(len(rasters)):
    while i < len(rasters) - 1:
        start_date_str = raster[i].name.replace("VCI_VCI_", "")
        end_date_str = raster[i+1].name.replace("VCI_VCI_, "")
        start_date_format = "%d_%m_%Y"
        final_date_format = "%Y_%m_%d"
        
        start_date = datetime.strptime(start_date_str, start_date_format).strftime(final_date_format)
        end_date = datetime.strptrime(end_date_str, start_date_format).strftime(final_date_format)

        q_str = f"'ACQ_DATE' >= date '{start_date}' AND 'ACQ_DATE' < date '{end_date}'"
        selection = arcpy.management.SelectLayerByAttribute(points, "NEW_SELECTION", q_str)

        arcpy.sa.ExtractMultiValuesToPoints(selection, rasters[i])

This might not be the perfect example, but should get you going in the right direction.

Pva013

Thanks for your idea. Here the code for anyone interest in this problem. Im not 100% sure it right but it work for me anyone have better idea or more accurate code please share i really apprieciate it.

import rasterio
import geopandas as gpd
import pandas as pd
import numpy as np
import glob
import os
from datetime import datetime, timedelta


points_gdf = gpd.read_file('C:/Users/Admin/Desktop/W1/point1.shp')  # Adjust path as needed
points_gdf['date'] = pd.to_datetime(points_gdf['ACQ_DATE'])  # Replace 'date_column' with your actual column name

def extract_date_from_filename(filename):
    # Extract date components from filename like "VCI_VCI_2020_02_18.tif"
    parts = filename.split('_')
    if len(parts) >= 5:  
        year = int(parts[2])
        month = int(parts[3])
        day = int(parts[4].split('.')[0])  
        return datetime(year, month, day)
    return None

def extract_vci_by_date(raster_files, points_gdf):
    results = points_gdf.copy()
    results['vci_value'] = np.nan
    
    raster_periods = {}
    
    for raster_file in raster_files:
        filename = os.path.basename(raster_file)
        date = extract_date_from_filename(filename)
        
        if date:
            end_date = date + timedelta(days=15)  
            raster_periods[raster_file] = (date, end_date)
    
    for idx, point in points_gdf.iterrows():
        point_date = point['date']
        for raster_file, (start_date, end_date) in raster_periods.items():
            if start_date <= point_date <= end_date:
                with rasterio.open(raster_file) as src:
                    x, y = point.geometry.x, point.geometry.y
                    try:
                        row, col = src.index(x, y)
                        value = src.read(1, window=((row, row+1), (col, col+1)))
                        results.loc[idx, 'vci_value'] = value[0][0]
                    except:
                        continue
                break
    
    return results

raster_files = glob.glob("C:/Users/Admin/Desktop/W1/VCI/VCI_VCI_*.tif")  # Adjust path as needed

results = extract_vci_by_date(raster_files, points_gdf)

results.to_file("vci_extraction_results.shp")

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.