卫星可以看到看不见的熔岩流和活跃野火，但如何做到？（Python）

原文towardsdatascience.com/satellites-can-see-invisible-lava-flows-and-active-wildires-but-how-python-371915464d1c目录 简介 Sentinel-2 (光谱波段)下载 Sentinel-2 图像⚙️处理 Sentinel-2 图像裁剪和重采样Sentinel-2 图像的可视化火山Sentinel-2 图像的可视化野火 结论 参考文献 简介如你所知我们的眼睛只能看到可见光区域蓝色、绿色和红色的波段。然而当光线击中物体并反射时它包含了其他光谱区域的信息如红外。红外光可以有效地穿透并穿过密集的气体如烟雾提供烟雾下的清晰视图。然而我们的眼睛无法看到红外区域的对象与一些动物如蛇可以看见它们视野中的一部分红外光不同。在过去几十年里红外光检测传感器的发展取得了显著进步。这些传感器已经用于实际应用中。我一直在寻找一个很好的例子来说明卫星如何检测红外区域的关键信息这些信息对我们肉眼来说是看不见的。上周我读到了关于冰岛火山自 2023 年 12 月以来第三次活跃的消息。这在我脑海中激发了一个想法那就是检查火山上卫星捕捉到的图像。我希望自己足够幸运能够找到一张清晰的火山烟雾柱卫星图像展示在可见光区域光线散射的同时红外区域穿透烟雾揭示熔岩流。我检查了几颗卫星猜猜看有一颗卫星的时机恰到好处火山在周四早上2 月 8 日爆发Sentinel-2 的飞越发生在 2 月 8 日中午左右。我认为这可以是一个完美的例子来展示卫星图像配备红外区域传感器如何帮助我们监测火山并检测到即使图像中有轻微烟雾也能发现熔岩流。为了探索在非常浓烟的情况下的红外波段的运用我决定可视化另一张 Sentinel-2 图像它捕捉到了 2020 年加利福尼亚州最大的野火之一。在这个故事中我们将下载两张 Sentinel-2 图像一张是在冰岛雷克雅内斯半岛上的火山上拍摄的另一张是在 2020 年加利福尼亚州 Creek 火灾中拍摄的。我们将使用 Python 通过可见光和红外区域的不同波段组合来可视化这两张图像。您将看到红外信息如何显示活跃的熔岩流和火点这些在红、绿、蓝RGB图像中是看不到的。如果您对此感兴趣请继续阅读 Sentinel-2 (光谱波段)Sentinel-2 任务由欧洲航天局ESA开发的 Sentinel-2A 和 Sentinel-2B 两颗卫星组成作为哥白尼计划的一部分。Sentinel-2 卫星配备了一种多光谱仪器可以捕捉 13 个光谱波段的数据。每个波段都有特定的波长范围允许进行广泛的地球观测应用。在本篇文章中我们将处理可见光区域的波段和红外区域的三个波段。以下是这些波段的列表波段 2蓝色496.6 纳米 波段 3绿色560.0 纳米 波段 4红色664.5 纳米波段 8近红外835.1 纳米 波段 11短波红外 11613.7 纳米 波段 12短波红外 22202.4 纳米波段 2蓝色、3绿色、4红色和 8近红外的空间分辨率为10 米波段 11 和 12短波红外为20 米这对于可视化火山和野火是足够的。下载 Sentinel-2 图像我已经为 Python 和 R 分别编写了两个下载 Sentinel-2 图像的教程。如果您想阅读完整的说明了解如何注册、设置凭据和下载图像请参阅以下帖子因为我将讨论在此帖中对图像下载代码进行必要调整使用 Google Colab 在 Python 中下载 Sentinel-2 图像更新于 2023 年 11 月如何在 R 中下载 Sentinel-2 图像更新于 2024 年 1 月根据那些脚本您需要修改感兴趣区域的坐标、图像日期并为两个示例火山和野火下载蓝色、绿色、红色、近红外和短波波段。每个案例都需要以下信息冰岛火山satelliteSENTINEL-2levelS2MSI1Caoi_pointPOINT(-22.411503 63.892295)start_date2024-02-07end_date2024-02-10Creek FiresatelliteSENTINEL-2levelS2MSI1Caoi_pointPOINT(-119.26 37.1914)start_date2020-09-07end_date2020-09-10一旦您获得了冰岛火山和 Creek 火灾可用图像的列表您将看到冰岛火山火山喷发发生在周四2 月 8 日早上 5:30 当地时间如图所示只有一张图像可用拍摄于 2 月 8 日 13:03 UTC12:30 当地时间正好在喷发后 7 小时。Creek Fire溪流火灾于 9 月 4 日开始于 12 月 24 日得到控制。列表中的两幅图像都是在 9 月 8 日拍摄的火灾开始后的第 4 天但第二幅图像的内容长度为零。因此我们将在这篇帖子中使用第一幅图像S2A_MSIL1C_20200908T183921_N0500_R070_T11SKB_20230309T124945。在运行下载 Sentinel-2 图像的这两个博客文章中的代码之前请确保在可见区域包含三个波段蓝色、绿色和红色和在红外区域包含三个波段近红外和短波红外波段。以下行对应这些波段f{product_name}/{root[0][0][12][0][0][1].text}.jp2# Bluef{product_name}/{root[0][0][12][0][0][2].text}.jp2# Greenf{product_name}/{root[0][0][12][0][0][3].text}.jp2# Redf{product_name}/{root[0][0][12][0][0][7].text}.jp2# Near-infraredf{product_name}/{root[0][0][12][0][0][11].text}.jp2# Shortwave infrared-1f{product_name}/{root[0][0][12][0][0][12].text}.jp2# Shortwave infrared-2如果你正确地遵循了步骤你应该在你的目录中为每个示例文件都有这些文件冰岛火山溪流火灾⚙️处理 Sentinel-2 图像剪切和降尺度在下载这些图像后我们需要剪切每个波段以适应火山和野火周围的感兴趣区域AOI。由于我们有每个事件的坐标我们可以使用此函数创建缓冲多边形火山为 3 公里野火为 10 公里importmathdefcalculate_new_coordinates(center_lat,center_lon,distance,bearing):# Earth radius in kilometersearth_radius6371.0# Convert coordinates to radianscenter_lat_radmath.radians(center_lat)center_lon_radmath.radians(center_lon)# Calculate new latitudenew_lat_radmath.asin(math.sin(center_lat_rad)*math.cos(distance/earth_radius)math.cos(center_lat_rad)*math.sin(distance/earth_radius)*math.cos(bearing))# Calculate new longitudenew_lon_radcenter_lon_radmath.atan2(math.sin(bearing)*math.sin(distance/earth_radius)*math.cos(center_lat_rad),math.cos(distance/earth_radius)-math.sin(center_lat_rad)*math.sin(new_lat_rad))# Convert back to degreesnew_latmath.degrees(new_lat_rad)new_lonmath.degrees(new_lon_rad)returnnew_lon,new_lat使用上面的函数进行冰岛火山AOI的处理# Center coordinatescenter_lat63.892295center_lon-22.411503# Buffer distance in kilometersbuffer_distance3# Calculate coordinates for the four cornersnorth_lon,north_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,0)south_lon,south_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,math.pi)east_lon,east_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,math.pi/2)west_lon,west_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,-math.pi/2)# Print the coordinates in the desired formatprint(f({west_lon:.4f}{north_lat:.4f},{east_lon:.4f}{north_lat:.4f},{east_lon:.4f}{south_lat:.4f},{west_lon:.4f}{south_lat:.4f},{west_lon:.4f}{north_lat:.4f}))# Output:(-22.472863.9193,-22.350263.9193,-22.350263.8653,-22.472863.8653,-22.472863.9193)使用上面的函数进行溪流火灾AOI的处理# Center coordinatescenter_lat37.19147center_lon-119.261175# Buffer distance in kilometersbuffer_distance10# Calculate coordinates for the four cornersnorth_lon,north_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,0)south_lon,south_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,math.pi)east_lon,east_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,math.pi/2)west_lon,west_latcalculate_new_coordinates(center_lat,center_lon,buffer_distance,-math.pi/2)# Print the coordinates in the desired formatprint(f({west_lon:.4f}{north_lat:.4f},{east_lon:.4f}{north_lat:.4f},{east_lon:.4f}{south_lat:.4f},{west_lon:.4f}{south_lat:.4f},{west_lon:.4f}{north_lat:.4f}))# Output:(-119.374137.2814,-119.148337.2814,-119.148337.1015,-119.374137.1015,-119.374137.2814)然后你可以使用以下坐标为 “aoi_polygon_wkt” 来剪切 jp2 文件正如我在本节 “堆叠转换为 Geotiff并剪切到 AOITOA” 中所描述的# Volcano: Polygon WKTaoi_polygon_wktPOLYGON ((-22.4728 63.9193, -22.3502 63.9193, -22.3502 63.8653, -22.4728 63.8653, -22.4728 63.9193))# Wildfire: Polygon WKTaoi_polygon_wktPOLYGON ((-119.3741 37.2814, -119.1483 37.2814, -119.1483 37.1015, -119.3741 37.1015, -119.3741 37.2814))使用这些多边形我们可以剪切 JP2 文件到我们的 AOI。然而在堆叠这些层进行可视化部分之前我们需要将短波红外波段从 20 米降尺度到 10 米以与其他波段蓝色、绿色、红色和近红外的大小兼容。这可以通过以下函数完成importrasteriofromscipy.ndimageimportzoomdefdownscale_raster(input_path,output_path,scale_factor):withrasterio.open(input_path)assrc:# Read the datadatasrc.read(1)# Calculate the new dimensionsnew_heightint(src.height/scale_factor)new_widthint(src.width/scale_factor)# Use scipys zoom function for resamplingresampled_datazoom(data,1/scale_factor,order3)# Update metadata for the new rastertransformsrc.transform*src.transform.scale((src.width/resampled_data.shape[1]),(src.height/resampled_data.shape[0]))new_profilesrc.profile new_profile.update({driver:JP2OpenJPEG,height:new_height,width:new_width,transform:transform})# Write the resampled raster to a new filewithrasterio.open(output_path,w,**new_profile)asdst:dst.write(resampled_data,1)然后我们可以应用此函数将短波波段从 20 米降尺度到 10 米以适应每个用例。使用上面的函数进行冰岛火山的处理# Usageinput_band_path_B11/content/T27VVL_20240208T130311_B11.jp2output_band_path_B11/content/T27VVL_20240208T130311_B11_resampled.jp2input_band_path_B12/content/T27VVL_20240208T130311_B12.jp2output_band_path_B12/content/T27VVL_20240208T130311_B12_resampled.jp2scale_factor1/2# 20m to 10mdownscale_raster(input_band_path_B11,output_band_path_B11,scale_factor)downscale_raster(input_band_path_B12,output_band_path_B12,scale_factor)使用上面的函数进行溪流火灾的处理# Usageinput_band_path_B11/content/T11SKB_20200908T183921_B11.jp2output_band_path_B11/content/T11SKB_20200908T183921_B11_resampled.jp2input_band_path_B12/content/T11SKB_20200908T183921_B12.jp2output_band_path_B12/content/T11SKB_20200908T183921_B12_resampled.jp2scale_factor1/2# 20m to 10mdownscale_raster(input_band_path_B11,output_band_path_B11,scale_factor)downscale_raster(input_band_path_B12,output_band_path_B12,scale_factor)在这些步骤之后你应该在你的目录中再有两个文件冰岛火山溪流火灾在所有波段都在同一维度的情况下我们可以根据本节 “堆叠转换为 Geotiff并剪切到 AOITOA” 中的内容继续前进生成堆叠文件。请参阅以下*文章*Sentinel-2 图像的可视化火山现在我们有两个堆叠的文件一个是火山事件文件另一个是野火事件文件我们可以使用不同的波段组合来绘制这些文件中的每一个。具体来说我们将创建三个图表一个仅基于可见波段红、绿、蓝另一个基于可见光和近红外波段绿、红和近红外的组合第三个专注于红外区域近红外和短波波段以了解如果我们省略红外区域的数据我们可能会错过哪些信息。要生成类似下面的图表你可以参考本*文章*中绘制 TOA 和地表反射率图像的真实颜色部分。冰岛火山这是在可见光区域记录的图像我们可以看到火山周围的熔岩黑色像素火山烟雾羽流以及一些非常小的红色区域显示着活跃的熔岩。如前所述光在可见波段容易被散射这就是为什么我们在这张图像中看到烟雾羽流为白色像素的原因。可见光范围内的光散射会遮挡物体使得观察烟雾下面的活跃熔岩变得具有挑战性。即使调整增益参数脚本中的增益参数控制亮度我们也只能看到向西流动的熔岩流在可见光谱中调整增益参数后你可以更清楚地看到向西流动的熔岩。然而我们仍然不知道烟雾羽流下面的情况。让我们再次绘制这张图像这次使用近红外信息。现在你可以看到当时有两股活跃的熔岩流一股向西流动在可见光范围内可以稍微检测到另一股向南流动通过近红外光揭示出来。向南流动的熔岩流靠近格里姆达维克镇该镇自去年 11 月那次喷发以来已被疏散。但问题是为什么即使在烟雾遮挡的情况下我们仍然能在红外区域看到熔岩我在引言中提到了这一点但背后的主要原因是红外光与可见光波段相比波长更长。这使得它能够穿过细小颗粒如密集气体和烟雾而不会发生散射与可见光波段不同。让我们继续前进再次绘制图像这次只使用红外波段具体是短波红外和近红外波段如你所见添加短波红外为图像增加了另一层信息。除了表示熔岩的亮度更高的像素外你还可以看到第一张和第二张图像中的黑色熔岩现在被分为两个区域红色和黑色。红色区域代表新燃烧的区域因为它在短波波段中反射更多可能含有热熔岩而其余部分显示的是非活跃熔岩。Sentinel-2 图像可视化野火与上一节类似我们将使用包括可见光波段红、蓝、绿、可见光和近红外波段绿、红和近红外以及红外波段短波和近红外的组合来绘制野火的 Sentinel-2 图像。让我们从可见光区域开始这是 Sentinel-2 在 Creek Fire 上捕获的可见光波段图像。如前所述光线在该区域容易散射这就是为什么唯一可见的是从燃烧区域上升到大气中的非常浓密的烟雾柱。然而正如火山示例所示红外区域的光线可以穿透烟雾揭示可见光区域隐藏的东西。为了评估浓烟中的情况让我们使用近红外和可见光波段绘制图像以揭示烟雾下面的情况没有太大区别如此所示与火山情况相比近红外类似于可见光波段已经散射对于揭示烟雾下面的情况没有帮助。由于 Sentinel-2 有两个比近红外波长更长的短波红外波段我们再次使用这两个短波红外波段和近红外波段绘制这张图像同时从可视化中移除任何可见光波段这不是很令人印象深刻吗短波红外波段能够穿透浓烟揭示出厚厚烟雾层下面的燃烧区域闪亮的金色像素这些区域仅凭可见光区域或甚至可见光和近红外区域的组合是无法看到的。定位这些燃烧区域对于火灾管理和野火监测可能是至关重要的。作为最后一步让我们将这些全部放在一起并使用此模板将这两个示例的图像并排绘制importrasterioimportnumpyasnpimportmatplotlib.pyplotasplt# Plot the stacked imagewithrasterio.open(stacked_TOA.tif)assrc:# Define band indicesblue_band1green_band2red_band3nir_band4swir1_band5swir2_band6# Read bandsredsrc.read(red_band)greensrc.read(green_band)bluesrc.read(blue_band)nirsrc.read(nir_band)swir1src.read(swir1_band)swir2src.read(swir2_band)# Apply gaingain2red_nnp.clip(red*gain/10000,0,1)green_nnp.clip(green*gain/10000,0,1)blue_nnp.clip(blue*gain/10000,0,1)nir_nnp.clip(nir*gain/10000,0,1)swir1_nnp.clip(swir1*gain/10000,0,1)swir2_nnp.clip(swir2*gain/10000,0,1)# Create different compositesrgb_compositenp.dstack((red_n,green_n,blue_n))nir_compositenp.dstack((nir_n,red_n,green_n))swir_compositenp.dstack((swir2_n,swir1_n,nir_n))# Plot the compositesplt.figure(figsize(24,8))plt.subplot(131)plt.title(Red, Green and Blue,fontsize18,fontweightbold)plt.imshow(rgb_composite)plt.axis(off)plt.subplot(132)plt.title(Near-infrared, Red, Green,fontsize18,fontweightbold)plt.imshow(nir_composite)plt.axis(off)plt.subplot(133)plt.title(Infrared Bands: Shortwave and Near-infrared,fontsize18,fontweightbold)plt.imshow(swir_composite)plt.axis(off)# Save the plotplt.savefig(composite_plot.png)plt.show()输出结果将是如在 Sentinel-2 部分所述这颗卫星有 13 个波段其中 9 个位于红外区域。我们只使用了红外区域的三波段来揭示烟雾下面的隐藏物体。您可以自由探索使用其他红外波段如红边波段波段 5、6、7 和 8A的视觉化并观察它如何影响结果。 结论图像之间的比较——一个在可见光区域捕获其他带有红外波段——揭示了在不同电磁谱中拥有眼睛和传感器的能力和重要性。额外的红外数据层使得能够清晰地识别和绘制出在 RGB 图像中不可见的活跃熔岩流和燃烧区域。与可见光相比红外波段的优越穿透力使得传感器即使在烟雾遮挡的情况下也能检测到活跃熔岩和活跃火灾等物体的存在。这种能力增强了我们评估这些事件潜在风险和及时做出灾害减轻决策的能力。 参考文献Copernicus Sentinel 数据 [2024] 用于 Sentinel 数据Copernicus 服务信息 [2024] 用于 Copernicus 服务信息。 在其他平台上与我联系获取更多精彩内容LinkedInResearchGateGithub以及Twitter。这里是通过以下链接可用的相关帖子从太空中观看风暴一个创建惊人视图的 Python 脚本如何在 R 中下载 Sentinel-2 图像更新至 2024 年 1 月使用 Google Colab 在 Python 中下载 Sentinel-2 图像更新至 2023 年 11 月使用 Python 和 Google Colab 可视化地理空间野火烟雾数据