用Python计算栅格数据的真实面积
在地理空间分析中，栅格数据的像素值通常代表某种属性，比如土地利用比例、植被覆盖率等。这些数据往往基于经纬度网格表示的比例值，而为了更直观地理解这些数据的空间意义，我们需要将这些比例值转化为实际面积（如平方米或公顷）。
对于高分辨率的大尺寸栅格文件（GeoTIFF），一次性加载整个文件会耗费大量内存，甚至导致处理失败。为了解决这个问题，我们可以通过分块处理 的方式，逐步计算像素的真实面积，确保即使是大文件也能顺利处理。
一. 应用场景
该方法尤其适用于以下场景：
1. 比例值栅格数据：
当栅格像素值表示某个区域的比例（例如森林覆盖率）时，通过面积计算，可以得出真实的覆盖面积。
2. 土地利用分析：
计算农田、草地、建筑用地等各类土地利用类型的真实面积。
3. 遥感影像处理：
在遥感影像中将像素值的属性映射到真实的地理空间。
二.代码实现
import sys
import os
import rasterio
import numpy as np
from glob import globdef haversine(coord1, coord2):"""Calculate the great circle distance in kilometers between two pointson the earth (specified in decimal degrees, with coordinates in the order of longitude, latitude).Arguments:coord1 -- tuple containing the longitude and latitude of the first point (lon1, lat1)coord2 -- tuple containing the longitude and latitude of the second point (lon2, lat2)Returns:distance in kilometers."""# Extract longitude and latitude, then convert from decimal degrees to radianslon1, lat1 = np.radians(coord1)lon2, lat2 = np.radians(coord2)# Haversine formuladlat = abs(lat2 - lat1)dlon = abs(lon2 - lon1)a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2c = 2 * np.arcsin(np.sqrt(a))r = 6371 # Radius of earth in kilometers. Use 3956 for milesreturn c * r# Read geotiff file in data folder
# Calculate area for a single GeoTIFF file
def calculate_area_for_tif(input_tif, chunk_size=1024):# Check if input file existsif not os.path.isfile(input_tif):raise FileNotFoundError(f"Input file '{input_tif}' does not exist.")# Construct output file pathoutput_path = input_tif.split('.')[0] + '_area.tif'# Skip if the area raster already existsif os.path.exists(output_path):print(f"Area of {input_tif} already calculated. Skipping...")return# Open the source rasterwith rasterio.open(input_tif) as src:# Check if the raster is in a geographic coordinate systemif not src.crs.is_geographic:raise ValueError(f"The raster {input_tif} is not in a geographic coordinate system!")# Get raster metadata and sizemeta = src.meta.copy()width, height = src.width, src.heighttransform = src.transform# Update metadata for outputmeta.update(compress='lzw', dtype=rasterio.float32, count=1)# Calculate the total number of chunkstotal_chunks = ((height + chunk_size - 1) // chunk_size) * ((width + chunk_size - 1) // chunk_size)chunk_count = 0 # Track processed chunks# Create the output raster filewith rasterio.open(output_path, 'w', **meta) as dst:print(f"Processing {input_tif} in chunks...")for row_start in range(0, height, chunk_size):row_end = min(row_start + chunk_size, height)for col_start in range(0, width, chunk_size):col_end = min(col_start + chunk_size, width)# Read a chunk of the source rasterwindow = rasterio.windows.Window(col_start, row_start, col_end - col_start, row_end - row_start)chunk_transform = src.window_transform(window)rows, cols = np.meshgrid(np.arange(row_start, row_end),np.arange(col_start, col_end),indexing='ij')# Apply geotransform to get geographical coordinatesx_coords, y_coords = chunk_transform * (cols, rows)x_coords_right, y_coords_right = chunk_transform * (cols + 1, rows)x_coords_bottom, y_coords_bottom = chunk_transform * (cols, rows + 1)xy_coords = np.stack((x_coords, y_coords), axis=0)xy_coords_right = np.stack((x_coords_right, y_coords_right), axis=0)xy_coords_bottom = np.stack((x_coords_bottom, y_coords_bottom), axis=0)# Calculate the area for this chunklength_right = haversine(xy_coords, xy_coords_right)length_bottom = haversine(xy_coords, xy_coords_bottom)area_chunk = length_right * length_bottom# Write the chunk to the output rasterdst.write(area_chunk.astype(np.float32), 1, window=window)# Update progresschunk_count += 1progress = (chunk_count / total_chunks) * 100print(f"Progress: {progress:.2f}% ({chunk_count}/{total_chunks} chunks)")print(f"\nArea of {input_tif} calculated and saved to {output_path}")if __name__ == '__main__':# Specify the input file pathinput_tif = 'path/to/your/input_file.tif'# Run the functioncalculate_area_for_tif(input_tif, chunk_size=1024)print('\nArea calculation complete!')
! 本代码修改自迪肯大学王金柱博士的代码
