import rasterio
from rasterio.windows import Window
import geopandas as gpd
import numpy as np
from joblib import Parallel, delayed
defprocess_block(data, window, transform, crs, target_crs, unit_factor):
"""
Process a single raster block for area calculation.
Parameters:
data (numpy array): Raster data block.
window (rasterio.windows.Window): Window corresponding to the data block.
transform (Affine): Affine transformation matrix for the data block.
crs (CRS): Original CRS of the raster.
target_crs (CRS): Target projection CRS for area calculation.
unit_factor (float): Conversion factor for the desired area unit.
Returns:
numpy array: Area calculation results for the block.
""" rows, cols = data.shape
polygons = []
for row inrange(rows):
for col inrange(cols):
x_min, y_min = transform * (col, row)
x_max, y_max = transform * (col + 1, row + 1)
polygon = {
'properties': {},
'geometry': {
'type': 'Polygon',
'coordinates': [[
[x_min, y_min],
[x_min, y_max],
[x_max, y_max],
[x_max, y_min],
[x_min, y_min]
]]
}
}
polygons.append(polygon)
gdf = gpd.GeoDataFrame.from_features(polygons, crs=crs)
ifstr(crs) != target_crs:
gdf = gdf.to_crs(target_crs)
gdf['area'] = gdf['geometry'].area * unit_factor
return np.array(gdf['area']).reshape(rows, cols)
defcalculate_pixel_area(input_tif, target_crs='EPSG:3577', area_unit='m', block_size=None, num_jobs=None):
"""
Calculate the area for each raster pixel and save the output as a raster file.
Parameters:
input_tif (str): Input raster file path.
target_crs (str): Target projection CRS for area calculation (default: 'EPSG:3577').
area_unit (str): Desired area unit ('m' for square meters, 'hm' for hectares, 'km' for square kilometers).
block_size (int or None): Block size for processing (e.g., 512). If None, processes the entire raster at once.
num_jobs (int or None): Number of parallel jobs for processing. Effective only when block_size is specified.
""" output_tif = input_tif.split('.')[0] + '_area.tif' unit_factors = {
'm': 1, # Square meters'hm': 1 / 10000, # Hectares'km': 1 / 1e6# Square kilometers }
if area_unit notin unit_factors:
raise ValueError("area_unit must be one of 'm', 'hm', or 'km'.")
with rasterio.open(input_tif) as src:
transform = src.transform
crs = src.crs
nodata = src.nodata if src.nodata isnotNoneelse0 dtype = 'float32'if block_size isNone:
# Process the entire raster without blocking data = src.read(1)
area_array = process_block(data, None, transform, crs, target_crs, unit_factors[area_unit])
else:
# Process in blocks height, width = src.height, src.width
windows = []
for row_off inrange(0, height, block_size):
for col_off inrange(0, width, block_size):
row_count = min(block_size, height - row_off)
col_count = min(block_size, width - col_off)
window = Window(col_off, row_off, col_count, row_count)
windows.append(window)
if num_jobs isNone:
# Single-threaded block processing results = [
process_block(
src.read(1, window=window),
window,
src.window_transform(window),
crs,
target_crs,
unit_factors[area_unit]
) for window in windows
]
else:
# Multi-threaded block processing results = Parallel(n_jobs=num_jobs)(
delayed(process_block)(
src.read(1, window=window),
window,
src.window_transform(window),
crs,
target_crs,
unit_factors[area_unit]
) for window in windows
)
# Combine block results into a full array area_array = np.zeros((height, width), dtype=dtype)
for result, window inzip(results, windows):
row_off, col_off = window.row_off, window.col_off
row_count, col_count = result.shape
area_array[row_off:row_off + row_count, col_off:col_off + col_count] = result
# Save the area raster filewith rasterio.open(
output_tif,
'w',
driver='GTiff',
height=src.height,
width=src.width,
count=1,
dtype=dtype,
crs=src.crs, # Preserve the original CRS transform=transform,
nodata=nodata
) as dst:
dst.write(area_array, 1)
print(f"Area raster file saved to: {output_tif}, unit: {area_unit}²")
if __name__ == '__main__':
# Input raster file path input_tif = 'rasters/Acrobates_pygmaeus.tif'# Call the function calculate_pixel_area(
input_tif=input_tif,
target_crs='EPSG:3577',
area_unit='km',
block_size=512, # Block size in pixels. None for no blocking. num_jobs=4# Number of parallel jobs, effective only when block_size is specified. )
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库