import os
import numpy as np
import rasterio
from rasterio.warp import Resampling
from rasterio.enums import Resampling as ResamplingEnum
from rasterio.merge import merge
from pathlib import Path
from skimage.measure import regionprops
from scipy.ndimage import label
import matplotlib.pyplot as plt
from matplotlib import colors
class Sentinel2AtmosphericCorrector:
"""Sentinel-2 L1C 大气校正处理器"""
def __init__(self, safe_path):
self.safe_path = Path(safe_path)
self.output_dir = self.safe_path.parent / f"{self.safe_path.stem}_CORRECTED"
os.makedirs(self.output_dir, exist_ok=True)
# 波段分辨率映射
self.band_resolutions = {
'B01': 60, 'B02': 10, 'B03': 10, 'B04': 10,
'B05': 20, 'B06': 20, 'B07': 20, 'B08': 10,
'B8A': 20, 'B09': 60, 'B10': 60, 'B11': 20, 'B12': 20
}
# 存储波段数据
self.band_data = {}
self.band_meta = {}
def find_band_path(self, band):
"""查找指定波段的文件路径"""
pattern = f"**/*{band}_{self.band_resolutions[band]}m.jp2"
return next(self.safe_path.glob(pattern))
def load_band(self, band):
"""加载指定波段数据"""
if band in self.band_data:
return self.band_data[band], self.band_meta[band]
band_path = self.find_band_path(band)
with rasterio.open(band_path) as src:
data = src.read(1).astype(np.float32)
meta = src.meta.copy()
self.band_data[band] = data / 10000 # 转换为反射率
self.band_meta[band] = meta
return self.band_data[band], self.band_meta[band]
def resample_band(self, target_band, reference_band):
"""将波段重采样到参考波段分辨率"""
target_data, _ = self.load_band(target_band)
ref_data, ref_meta = self.load_band(reference_band)
# 如果需要重采样
if self.band_resolutions[target_band] != self.band_resolutions[reference_band]:
scale_factor = self.band_resolutions[target_band] / self.band_resolutions[reference_band]
new_shape = (int(target_data.shape[0] * scale_factor),
int(target_data.shape[1] * scale_factor))
# 使用rasterio进行高质量重采样
with rasterio.MemoryFile() as memfile:
with memfile.open(
driver='GTiff',
height=target_data.shape[0],
width=target_data.shape[1],
count=1,
dtype=np.float32
) as temp_ds:
temp_ds.write(target_data, 1)
with memfile.open() as temp_ds:
resampled, _ = reproject(
source=rasterio.band(temp_ds, 1),
destination=np.empty(new_shape, dtype=np.float32),
src_transform=temp_ds.transform,
src_crs=temp_ds.crs,
dst_transform=ref_meta['transform'],
dst_crs=ref_meta['crs'],
resampling=Resampling.bilinear
)
return resampled
return target_data
def auto_detect_targets(self):
"""自动检测暗目标和亮目标区域"""
# 1. 加载关键波段
b08, _ = self.load_band('B08') # 近红外
b11, _ = self.load_band('B11') # 短波红外
# 将B11重采样到10m分辨率
b11 = self.resample_band('B11', 'B08')
# 2. 检测暗目标（水体）
water_mask = (b08 > 0.01) & (b08 < 0.1)
labeled_water, num_features = label(water_mask)
regions = regionprops(labeled_water)
if not regions:
raise ValueError("未检测到适合的水体区域作为暗目标")
# 选择最大的三个水体区域
large_waters = sorted(regions, key=lambda r: r.area, reverse=True)[:3]
centers = np.array([r.centroid for r in large_waters])
avg_center = np.mean(centers, axis=0).astype(int)
# 创建暗目标区域（100x100像素）
roi_size = 100
dark_coords = (
max(0, avg_center[1] - roi_size//2),
max(0, avg_center[0] - roi_size//2),
min(b08.shape[1], avg_center[1] + roi_size//2),
min(b08.shape[0], avg_center[0] + roi_size//2)
)
# 3. 检测亮目标（云层或沙地）
cloud_mask = (b11 > 0.8)
labeled_clouds, _ = label(cloud_mask)
cloud_regions = regionprops(labeled_clouds)
if cloud_regions:
# 选择最大的三个云层区域
large_clouds = sorted(cloud_regions, key=lambda r: r.area, reverse=True)[:3]
centers = np.array([r.centroid for r in large_clouds])
avg_center = np.mean(centers, axis=0).astype(int)
else:
# 若无云层则使用沙地区域
sand_mask = (b08 > 0.3) & (b08 < 0.6)
labeled_sand, _ = label(sand_mask)
sand_regions = regionprops(labeled_sand)
if not sand_regions:
raise ValueError("未检测到合适的亮目标区域")
large_sands = sorted(sand_regions, key=lambda r: r.area, reverse=True)[:3]
centers = np.array([r.centroid for r in large_sands])
avg_center = np.mean(centers, axis=0).astype(int)
# 创建亮目标区域
bright_coords = (
max(0, avg_center[1] - roi_size//2),
max(0, avg_center[0] - roi_size//2),
min(b08.shape[1], avg_center[1] + roi_size//2),
min(b08.shape[0], avg_center[0] + roi_size//2)
)
return dark_coords, bright_coords
def extract_roi(self, data, coords):
"""提取感兴趣区域数据"""
xmin, ymin, xmax, ymax = coords
return data[ymin:ymax, xmin:xmax]
def calculate_correction_params(self, band):
"""计算大气校正参数"""
# 获取波段数据
band_data, _ = self.load_band(band)
# 获取目标区域
dark_roi = self.extract_roi(band_data, self.dark_coords)
bright_roi = self.extract_roi(band_data, self.bright_coords)
# 计算目标区域平均值
dark_mean = np.nanmean(dark_roi)
bright_mean = np.nanmean(bright_roi)
# 假设反射率值
dark_refl = 0.01 # 水体典型反射率
bright_refl = 0.85 # 云层或沙地典型反射率
# 计算线性校正参数
slope = (bright_refl - dark_refl) / (bright_mean - dark_mean)
intercept = dark_refl - slope * dark_mean
return slope, intercept
def apply_correction(self):
"""执行大气校正并保存结果"""
# 1. 自动检测目标区域
self.dark_coords, self.bright_coords = self.auto_detect_targets()
print(f"📊 自动检测结果:")
print(f" - 暗目标区域坐标: {self.dark_coords}")
print(f" - 亮目标区域坐标: {self.bright_coords}")
# 2. 可视化目标区域
self.visualize_targets()
# 3. 处理所有波段
processed_bands = {}
for band in self.band_resolutions.keys():
try:
# 加载波段数据
band_data, meta = self.load_band(band)
# 计算校正参数
slope, intercept = self.calculate_correction_params(band)
# 应用校正
corrected = slope * band_data + intercept
corrected = np.clip(corrected, 0, 1) # 限制在0-1范围
# 保存结果
output_path = self.output_dir / f"{band}_CORRECTED.tif"
meta.update({
'dtype': 'float32',
'driver': 'GTiff',
'count': 1
})
with rasterio.open(output_path, 'w', **meta) as dst:
dst.write(corrected.astype(np.float32), 1)
print(f"✅ 波段 {band} 校正完成 -> {output_path}")
processed_bands[band] = corrected
except Exception as e:
print(f"⚠️ 波段 {band} 处理失败: {str(e)}")
# 4. 生成RGB预览图
self.generate_rgb_preview(processed_bands)
return processed_bands
def visualize_targets(self):
"""可视化目标区域位置"""
b08, meta = self.load_band('B08')
# 创建可视化图像
fig, ax = plt.subplots(figsize=(12, 12))
ax.imshow(b08, cmap='gray', vmin=0, vmax=0.3)
# 绘制暗目标区域
rect = plt.Rectangle(
(self.dark_coords[0], self.dark_coords[1]),
self.dark_coords[2] - self.dark_coords[0],
self.dark_coords[3] - self.dark_coords[1],
linewidth=2, edgecolor='blue', facecolor='none'
)
ax.add_patch(rect)
ax.text(self.dark_coords[0], self.dark_coords[1]-20,
'Dark Target', color='blue', fontsize=12,
bbox=dict(facecolor='white', alpha=0.7))
# 绘制亮目标区域
rect = plt.Rectangle(
(self.bright_coords[0], self.bright_coords[1]),
self.bright_coords[2] - self.bright_coords[0],
self.bright_coords[3] - self.bright_coords[1],
linewidth=2, edgecolor='red', facecolor='none'
)
ax.add_patch(rect)
ax.text(self.bright_coords[0], self.bright_coords[1]-20,
'Bright Target', color='red', fontsize=12,
bbox=dict(facecolor='white', alpha=0.7))
plt.title(f"目标区域检测结果 - {self.safe_path.stem}", fontsize=14)
plt.axis('off')
output_path = self.output_dir / "target_detection.png"
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"🌄 目标区域可视化已保存至: {output_path}")
def generate_rgb_preview(self, bands):
"""生成RGB预览图像"""
try:
# 获取RGB波段
red = bands['B04']
green = bands['B03']
blue = bands['B02']
# 创建RGB合成
rgb = np.stack([red, green, blue], axis=-1)
# 应用拉伸增强
p2, p98 = np.percentile(rgb, (2, 98))
rgb_stretched = np.clip((rgb - p2) / (p98 - p2), 0, 1)
fig, ax = plt.subplots(figsize=(15, 15))
ax.imshow(rgb_stretched)
plt.title(f"大气校正结果预览 - {self.safe_path.stem}", fontsize=16)
plt.axis('off')
output_path = self.output_dir / "corrected_rgb_preview.png"
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"🖼️ RGB预览图已保存至: {output_path}")
except Exception as e:
print(f"⚠️ 无法生成RGB预览: {str(e)}")
# 使用示例
if __name__ == "__main__":
# 指定影像路径
safe_path = r"D:\大三上\经验线性校正\S2B_MSIL1C_20181009T110939_N0500_R137_T29TQG_20230714T221126.SAFE"
# 创建校正器并执行校正
corrector = Sentinel2AtmosphericCorrector(safe_path)
corrected_bands = corrector.apply_correction()
print(f"🎉 大气校正完成! 结果保存在: {corrector.output_dir}")
代码如上，报错如下A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.3.4 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "<frozen runpy>", line 198, in _run_module_as_main
File "<frozen runpy>", line 88, in _run_code
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel_launcher.py", line 18, in <module>
app.launch_new_instance()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\traitlets\config\application.py", line 1075, in launch_instance
app.start()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\kernelapp.py", line 739, in start
self.io_loop.start()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\tornado\platform\asyncio.py", line 195, in start
self.asyncio_loop.run_forever()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\asyncio\base_events.py", line 608, in run_forever
self._run_once()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\asyncio\base_events.py", line 1936, in _run_once
handle._run()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\asyncio\events.py", line 84, in _run
self._context.run(self._callback, *self._args)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\kernelbase.py", line 545, in dispatch_queue
await self.process_one()
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\kernelbase.py", line 534, in process_one
await dispatch(*args)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\kernelbase.py", line 437, in dispatch_shell
await result
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\ipkernel.py", line 362, in execute_request
await super().execute_request(stream, ident, parent)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\kernelbase.py", line 778, in execute_request
reply_content = await reply_content
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\ipkernel.py", line 449, in do_execute
res = shell.run_cell(
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\ipykernel\zmqshell.py", line 549, in run_cell
return super().run_cell(*args, **kwargs)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\interactiveshell.py", line 3048, in run_cell
result = self._run_cell(
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\interactiveshell.py", line 3103, in _run_cell
result = runner(coro)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\async_helpers.py", line 129, in _pseudo_sync_runner
coro.send(None)
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\interactiveshell.py", line 3308, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\interactiveshell.py", line 3490, in run_ast_nodes
if await self.run_code(code, result, async_=asy):
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\IPython\core\interactiveshell.py", line 3550, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "C:\Users\ywy\AppData\Local\Temp\ipykernel_5116\349592770.py", line 10, in <module>
import matplotlib.pyplot as plt
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\__init__.py", line 113, in <module>
from . import _api, _version, cbook, _docstring, rcsetup
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\rcsetup.py", line 27, in <module>
from matplotlib.colors import Colormap, is_color_like
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\colors.py", line 56, in <module>
from matplotlib import _api, _cm, cbook, scale
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\scale.py", line 22, in <module>
from matplotlib.ticker import (
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\ticker.py", line 138, in <module>
from matplotlib import transforms as mtransforms
File "D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\transforms.py", line 49, in <module>
from matplotlib._path import (
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
AttributeError: _ARRAY_API not found
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[16], line 10
8 from skimage.measure import regionprops
9 from scipy.ndimage import label
---> 10 import matplotlib.pyplot as plt
11 from matplotlib import colors
13 class Sentinel2AtmosphericCorrector:
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\__init__.py:113
109 from packaging.version import parse as parse_version
111 # cbook must import matplotlib only within function
112 # definitions, so it is safe to import from it here.
--> 113 from . import _api, _version, cbook, _docstring, rcsetup
114 from matplotlib.cbook import sanitize_sequence
115 from matplotlib._api import MatplotlibDeprecationWarning
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\rcsetup.py:27
25 from matplotlib import _api, cbook
26 from matplotlib.cbook import ls_mapper
---> 27 from matplotlib.colors import Colormap, is_color_like
28 from matplotlib._fontconfig_pattern import parse_fontconfig_pattern
29 from matplotlib._enums import JoinStyle, CapStyle
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\colors.py:56
54 import matplotlib as mpl
55 import numpy as np
---> 56 from matplotlib import _api, _cm, cbook, scale
57 from ._color_data import BASE_COLORS, TABLEAU_COLORS, CSS4_COLORS, XKCD_COLORS
60 class _ColorMapping(dict):
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\scale.py:22
20 import matplotlib as mpl
21 from matplotlib import _api, _docstring
---> 22 from matplotlib.ticker import (
23 NullFormatter, ScalarFormatter, LogFormatterSciNotation, LogitFormatter,
24 NullLocator, LogLocator, AutoLocator, AutoMinorLocator,
25 SymmetricalLogLocator, AsinhLocator, LogitLocator)
26 from matplotlib.transforms import Transform, IdentityTransform
29 class ScaleBase:
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\ticker.py:138
136 import matplotlib as mpl
137 from matplotlib import _api, cbook
--> 138 from matplotlib import transforms as mtransforms
140 _log = logging.getLogger(__name__)
142 __all__ = ('TickHelper', 'Formatter', 'FixedFormatter',
143 'NullFormatter', 'FuncFormatter', 'FormatStrFormatter',
144 'StrMethodFormatter', 'ScalarFormatter', 'LogFormatter',
(...)
150 'MultipleLocator', 'MaxNLocator', 'AutoMinorLocator',
151 'SymmetricalLogLocator', 'AsinhLocator', 'LogitLocator')
File D:\ArcGISpro3.4.3\bin\Python\envs\arcgispro-py3\Lib\site-packages\matplotlib\transforms.py:49
46 from numpy.linalg import inv
48 from matplotlib import _api
---> 49 from matplotlib._path import (
50 affine_transform, count_bboxes_overlapping_bbox, update_path_extents)
51 from .path import Path
53 DEBUG = False
ImportError: numpy.core.multiarray failed to import
怎么改
本文介绍了解决在使用skimage的measure模块时遇到的label报错问题，该问题通常由于skimage版本过低导致。文章提供了一个通用的解决方案，即通过指定清华的源来升级skimage，确保代码能够正常运行。
扫一扫
726
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
