基于非线性函数的采样点周边多边形面积加权和与栅格均值计算

Got it, let's walk through how to solve this problem step by step. You've got 100 sampling points, plus polygon and raster environmental background data, and need to compute distance-decay weighted metrics for each point—here's a practical breakdown tailored to both data types:

The key here is using your nonlinear distance-decay function to assign higher weights to environmental features closer to each sampling point, and lower weights to those farther away. For context, I'll assume your function follows a common form like w(r) = exp(-kr) (exponential decay) or w(r) = 1/(1+kr²) (inverse squared decay), where r is the distance from the sampling point to the feature.

1. Polygon Data: Calculate Area-Weighted Sum #

This method accounts for both the distance of each polygon segment and its size when computing the weighted sum of your environmental variable.

• Step 1: Define Influence Zones
  First, set a maximum distance threshold (e.g., 5km) based on your research needs—only polygons within this buffer around a sampling point will be considered. Use spatial tools to generate buffers for each point, then clip your environmental polygons to these buffers to get valid, overlapping segments.

  • For Python users: Use geopandas.GeoSeries.buffer() and geopandas.clip()
  • For GIS desktop users: ArcGIS Buffer + Clip tools, QGIS Buffer + Clip

• Step 2: Compute Distance Weights
  For each clipped polygon segment, calculate the distance r from the segment to the sampling point. A common shortcut is using the segment's centroid distance, but you can use minimum distance (closest edge) for higher precision. Plug r into your nonlinear decay function to get the weight w(r) for that segment.

• Step 3: Calculate Weighted Sum
  For your target environmental variable, compute the product of the variable value, segment area, and weight for each segment. Sum these products to get the final area-weighted sum for the sampling point.

Python Code Example (GeoPandas) #

import geopandas as gpd
import numpy as np

# Load your data
gdf_points = gpd.read_file("sampling_points.shp")
gdf_polygons = gpd.read_file("environmental_polygons.shp")

# Define your distance-decay function (replace with your actual formula)
def distance_decay(r, k=0.001):
    # Example: Exponential decay, k adjusts the rate of decay
    return np.exp(-k * r)

# Process each sampling point
for idx, point in gdf_points.iterrows():
    # Create 1km buffer around the point (adjust distance as needed)
    buffer = point.geometry.buffer(1000)
    # Clip polygons to the buffer
    clipped_segments = gpd.clip(gdf_polygons, buffer)
    
    if clipped_segments.empty:
        gdf_points.loc[idx, "polygon_weighted_sum"] = 0
        continue
    
    # Calculate centroid distance from each segment to the sampling point
    clipped_segments["distance"] = clipped_segments.geometry.centroid.distance(point.geometry)
    # Compute weights
    clipped_segments["weight"] = clipped_segments["distance"].apply(distance_decay)
    # Calculate segment area (in square meters)
    clipped_segments["area"] = clipped_segments.geometry.area
    
    # Compute weighted sum: sum(variable * area * weight)
    weighted_sum = (clipped_segments["env_variable"] * clipped_segments["area"] * clipped_segments["weight"]).sum()
    gdf_points.loc[idx, "polygon_weighted_sum"] = weighted_sum

# Save results
gdf_points.to_file("weighted_points_polygon.shp")

2. Raster Data: Calculate Weighted Mean #

For raster data, we'll weight each pixel based on its distance to the sampling point, then compute a weighted average of the pixel values.

• Step 1: Extract Relevant Pixels
  Use the same buffer threshold as before to clip your raster to the area around each sampling point. This isolates only the pixels that contribute to the weighted mean.

  • For Python users: rasterio.mask.mask()
  • For GIS desktop users: ArcGIS Extract by Mask, QGIS Clip Raster by Mask Layer

• Step 2: Compute Pixel Weights
  For each clipped pixel, calculate the distance r from its center to the sampling point. Apply your decay function to get w(r), making sure to ignore NoData pixels.

• Step 3: Compute Weighted Mean
  The weighted mean is the sum of (pixel value × weight) divided by the sum of all weights. This accounts for both the pixel's value and its relative importance based on distance.

Python Code Example (Rasterio) #

import rasterio
from rasterio.mask import mask
import numpy as np
from shapely.geometry import Point

# Load raster and sampling points
with rasterio.open("environmental_raster.tif") as src:
    # Assume sampling points are stored as a list of (x,y) tuples
    sampling_points = [(x, y) for x, y in zip(gdf_points.x, gdf_points.y)]
    weighted_means = []
    
    for x, y in sampling_points:
        # Create buffer geometry
        point_geom = Point(x, y)
        buffer_geom = point_geom.buffer(1000)
        
        # Clip raster to buffer
        out_image, out_transform = mask(src, [buffer_geom], crop=True, nodata=src.nodata)
        # Flatten raster to 1D array, filter out NoData
        pixel_values = out_image[0][out_image[0] != src.nodata]
        
        if len(pixel_values) == 0:
            weighted_means.append(0)
            continue
        
        # Calculate coordinates of each pixel center
        rows, cols = np.where(out_image[0] != src.nodata)
        xs, ys = rasterio.transform.xy(out_transform, rows, cols)
        # Compute distance from each pixel to sampling point
        distances = np.sqrt((np.array(xs) - x)**2 + (np.array(ys) - y)**2)
        # Calculate weights using your decay function
        weights = distance_decay(distances)
        
        # Compute weighted mean
        weighted_mean = np.sum(pixel_values * weights) / np.sum(weights)
        weighted_means.append(weighted_mean)

# Add results to sampling points DataFrame
gdf_points["raster_weighted_mean"] = weighted_means
gdf_points.to_file("weighted_points_raster.shp")

Key Tips for Success #

• Calibrate Decay Parameters: The k value (or other parameters in your nonlinear function) should be tuned to your study context—check existing literature in your field for typical values, or use field data to validate.
• Optimize Performance: For 100 points, looping is manageable, but for larger datasets, consider batch processing tools or spatial indexing (e.g., geopandas.sjoin() to pre-filter polygons near points).
• Distance Precision: If your research requires high accuracy, use minimum distance (polygon edge to point) instead of centroid distance, or pixel edge distance instead of center distance.

内容的提问来源于stack exchange，提问作者Nick_89