如何给网站配色,国内设计大神网站,企业怎么做app网址,沈阳建设工程城乡建设厅将 Segment Anything Model Version 2 应用于卫星图像以检测和导出农业地区田地边界的分步教程
#x1f31f; 简介
手动绘制田地边界是最耗时的任务之一#xff0c;其准确性取决于绘制者的表现。然而#xff0c;精确的边界检测在很多领域都有应用。例如#xff0c;假设您…将 Segment Anything Model Version 2 应用于卫星图像以检测和导出农业地区田地边界的分步教程 简介
手动绘制田地边界是最耗时的任务之一其准确性取决于绘制者的表现。然而精确的边界检测在很多领域都有应用。例如假设您想训练一种机器学习算法分析卫星图像中的植被指数与农场作物产量之间的关系。您需要的第一个输入是农场的形状文件这通常需要手动绘制。绘制一个形状文件可能只需要几分钟但如果您需要为 1000 个农场绘制边界呢这时这个过程就变得非常耗时而自动提取边界的技术就变得非常有价值可以节省数小时的工作时间。
在本教程中我将演示如何使用由吴秋生博士基于第一版和第二版 分段任何模型SAM开发的 segment-anything-py 和 segment-geospatial Python 软件包。所有代码都是在 Google Colab 中编写和测试的任何人都可以轻松复制这些步骤。如果您对此感兴趣请继续阅读 设置 Google Colab
所有代码都将使用 Python 编写并在 Google Colab 平台上进行测试因此您无需安装各种软件和编译器即可按照步骤进行操作。由于运行 SAM 需要 GPU因此请确保将运行时更改为 TPUv4方法是点击 运行时 选项卡选择 更改运行时类型然后选择 TPUv4。此外还需要使用 pip 命令安装以下软件包
pip install pandas rasterio
️ 加载清晰的哨兵-2 图像
设置好 Google Colab 后我们需要一张农田的航空图像。我在本教程中使用了一张哨兵-2 图像但您也可以使用任何按顺序蓝、绿、红保存了红色、绿色和蓝色波段的卫星图像。
Downloading Sentinel-2 Imagery in Python with Google Colab (Updated Nov 2023)
并使用以下信息检索相同的图像图像信息 (S2B_MSIL2A_20240806T184919_N0511_R113_T10SFH)
url_dataspace https://catalogue.dataspace.copernicus.eu/odata/v1satellite SENTINEL-2
level S2MSI2Aaoi_point POINT(-121.707902 38.368628)cloud_cover 10start_date 2024-07-15
end_date 2024-08-10
start_date_full start_dateT00:00:00.000Z
end_date_full end_date T00:00:00.000Z 按照这些步骤操作后您的内容文件夹中就会出现 JP2 格式的三个单独色带红、绿、蓝如下图所示 在哨兵-2 图像上应用 SAM2
将 SAM2 应用于卫星图像相对简单但需要额外的步骤为模型准备图像。第一步是剪切下载的场景将重点放在我们感兴趣的区域AOI上因为完整的场景可能包括我们不感兴趣的区域如城区、海洋、湖泊、山脉或森林。此外Google Colab 的资源可能不足以处理整个场景。要创建一个较小的 AOI我们可以在农业区域内定义一个点并在该点的周围设置一个约 5 千米的缓冲区。
第二步是保存剪切后的图像并将波段排序为蓝、绿、红BGR因为算法希望采用这种顺序而不是通常的 RGB。最后将输出保存为 GeoTIFF 格式因为算法不接受 JP2 格式的文件。下面的代码在点周围定义了一个缓冲区根据边界框剪切红、绿、蓝三色带并以 BGR 顺序将输出保存为 GeoTIFF 格式
import rasterio
from rasterio.merge import merge
from rasterio.plot import show
from rasterio.mask import mask
from shapely.geometry import Point, box
from shapely.wkt import loads as load_wkt
import geopandas as gpd
from pyproj import CRS, Transformer
import numpy as np
import osdef clip_and_merge_jp2_files(blue_jp2, green_jp2, red_jp2, aoi_point_wkt, buffer_radius_km, output_tiff):# Parse the AOI point from WKTaoi_point load_wkt(aoi_point_wkt)# Open the JP2 fileswith rasterio.open(blue_jp2) as blue_src, \rasterio.open(green_jp2) as green_src, \rasterio.open(red_jp2) as red_src:# Get the CRS of the JP2 files jp2_crs blue_src.crs# Create a GeoDataFrame for the AOI point aoi_gdf gpd.GeoDataFrame({geometry: [aoi_point]}, crsEPSG:4326)# Reproject the AOI point to the JP2 CRS if aoi_gdf.crs ! jp2_crs:aoi_gdf aoi_gdf.to_crs(jp2_crs)# Create a buffer around the AOI point (in meters)buffer_radius buffer_radius_km * 1000 # Convert km to metersaoi_buffer aoi_gdf.geometry.buffer(buffer_radius).iloc[0]# Convert the buffer to a bounding boxminx, miny, maxx, maxy aoi_buffer.boundsbbox box(minx, miny, maxx, maxy)# Convert the bbox to a GeoDataFramebbox_gdf gpd.GeoDataFrame({geometry: [bbox]}, crsjp2_crs)# Clip each band using the bboxblue_clipped, blue_transform mask(blue_src, bbox_gdf.geometry, cropTrue)green_clipped, green_transform mask(green_src, bbox_gdf.geometry, cropTrue)red_clipped, red_transform mask(red_src, bbox_gdf.geometry, cropTrue)# Update the metadata meta blue_src.meta.copy()meta.update({driver: GTiff,height: blue_clipped.shape[1],width: blue_clipped.shape[2],transform: blue_transform,count: 3, # We have three bands: B, G, Rdtype: blue_clipped.dtype})# Merge the bands into a single arraymerged_bgr np.stack([blue_clipped[0], green_clipped[0], red_clipped[0]])# Save the merged BGR image as a GeoTIFFwith rasterio.open(output_tiff, w, **meta) as dst:dst.write(merged_bgr)print(fClipped and merged image saved as {output_tiff})blue_jp2 T10SFH_20240806T184919_B02_10m.jp2
green_jp2 T10SFH_20240806T184919_B03_10m.jp2
red_jp2 T10SFH_20240806T184919_B04_10m.jp2
buffer_radius_km 1.5
output_tiff BGR_20240806.tif
aoi_point POINT(-121.707902 38.368628) #AOI point (longitude, latitude)clip_and_merge_jp2_files(blue_jp2, green_jp2, red_jp2, aoi_point, buffer_radius_km, output_tiff)
运行代码后你应该能在内容文件夹中看到剪切后的图片
import matplotlib.pyplot as pltdef plot_tiff(tiff_file):# Open the tiff filewith rasterio.open(tiff_file) as src:b_band src.read(1) g_band src.read(2) r_band src.read(3) # Stack the bands into a single numpy arrayrgb np.dstack((r_band, g_band, b_band))# Normalize the bands to the range [0, 1] (for display)rgb rgb.astype(np.float32)rgb / np.max(rgb)# Plot the imageplt.imshow(rgb)plt.axis(off) # Hide the axisplt.show()plot_tiff(BGR_20240806.tif) 下载图像后下一步是剪切图像并将其保存为可接受的格式。我们需要更改图像格式因为算法需要 8 位无符号格式而剪切后的图像是浮点格式。下面的脚本转换了格式并以 8 位无符号格式保存图像
def convert_to_8bit(input_tiff, output_tiff):with rasterio.open(input_tiff) as src:blue src.read(1)green src.read(2)red src.read(3)# Normalize the float values to 0-255 and convert to 8-bit unsigned integersblue_8bit np.clip((blue - np.min(blue)) / (np.max(blue) - np.min(blue)) * 255, 0, 255).astype(np.uint8)green_8bit np.clip((green - np.min(green)) / (np.max(green) - np.min(green)) * 255, 0, 255).astype(np.uint8)red_8bit np.clip((red - np.min(red)) / (np.max(red) - np.min(red)) * 255, 0, 255).astype(np.uint8)# Define metadata profile src.profileprofile.update(dtyperasterio.uint8,count3,compresslzw)# Write the new 8-bit data to the output filewith rasterio.open(output_tiff, w, **profile) as dst:dst.write(blue_8bit, 1)dst.write(green_8bit, 2)dst.write(red_8bit, 3)input_tiff BGR_20240806.tif
output_tiff BGR_20240806_8bit.tif
convert_to_8bit(input_tiff, output_tiff) 第三步是将 UTM 坐标的图像保存为地理坐标经纬度。运行以下代码即可完成此操作
from rasterio.warp import calculate_default_transform, reproject, Resamplingdef convert_to_latlong(input_tiff, output_tiff):with rasterio.open(input_tiff) as src:transform, width, height calculate_default_transform(src.crs, EPSG:4326, src.width, src.height, *src.bounds)kwargs src.meta.copy()kwargs.update({crs: EPSG:4326,transform: transform,width: width,height: height})with rasterio.open(output_tiff, w, **kwargs) as dst:for i in range(1, src.count 1):reproject(sourcerasterio.band(src, i),destinationrasterio.band(dst, i),src_transformsrc.transform,src_crssrc.crs,dst_transformtransform,dst_crsEPSG:4326,resamplingResampling.nearest)input_tiff BGR_20240806.tif
output_tiff BGR_20240806_reproj.tif
convert_to_latlong(input_tiff, output_tiff)
最后一步取决于您想如何部署和使用 SAM 算法。有两种模式可供选择自动和手动。在自动模式下算法只需要我们导出的准备好的图像带地理坐标的 8 位无符号格式剪切 BGR 图像。在手动模式下您可以在每个对象上添加一个点这通常有助于算法生成更精确的结果并对用户点识别的对象进行分割。要在自动模式下运行算法可以跳过下面的章节直接跳到 自动模式下的 SAM。但是如果您还想使用手动模式请添加下面的脚本这样您就可以点击图像并以经纬度存储您的点。
from localtileserver import get_folium_tile_layer, TileClient,get_leaflet_tile_layer
import ipyleaflet
from shapely.geometry import Point
from ipyleaflet import Map, Marker, ImageOverlay
from ipywidgets import Output, VBox
from IPython.display import display
import matplotlib.pyplot as plt
from PIL import Imagegeotiff_path BGR_20240806_reproj.tif# Create a TileClient object
client TileClient(geotiff_path)# Create a TileLayer using the client
tiff_layer get_leaflet_tile_layer(client, nameGeoTIFF)# Get the bounds of the GeoTIFF
bounds client.bounds()
center ((bounds[0] bounds[1]) / 2, (bounds[2] bounds[3]) / 2)# Create an ipyleaflet map
m Map(centercenter, zoom14)# Add the TileLayer to the map
m.add_layer(tiff_layer)# Create a list to store the clicked points
clicked_points []# Create an output widget to capture map click events
output Output()# Function to handle clicks on the map
def handle_click(**kwargs):if type in kwargs and kwargs[type] click:latlon kwargs.get(coordinates)if latlon:lat, lon latlonclicked_points.append(Point(lon, lat))marker Marker(location(lat, lon))m.add_layer(marker)with output:print(fPoint added: {lat}, {lon})# Add the click handler to the map
m.on_interaction(handle_click)# Display the map and output widget
display(VBox([m, output]))运行代码后会出现一张交互式地图您可以点击地图。每次点击后这些点都会用蓝色标记标出如下图所示要查看您在地图上所选点的坐标只需运行以下代码即可
clicked_points
[POINT (-121.709 38.371),POINT (-121.716 38.371),POINT (-121.717 38.37),POINT (-121.717 38.368),POINT (-121.717 38.366),POINT (-121.709 38.366),POINT (-121.709 38.369),POINT (-121.7 38.371),POINT (-121.701 38.369),POINT (-121.7 38.367),POINT (-121.697 38.375),POINT (-121.715 38.377),POINT (-121.718 38.379),POINT (-121.72 38.363),POINT (-121.699 38.362)]
您还可以通过使用
# Function to export the points to a GeoPackage
def export_to_gpkg(points, output_path):Export points to a GeoPackage.gdf gpd.GeoDataFrame(geometrypoints, crsEPSG:4326)gdf.to_file(output_path, driverGPKG)output_gpkg_path output.gpkg
export_to_gpkg(clicked_points, output_gpkg_path)
自动模式的 SAM
如前所述如果输入图像的格式符合 SAM 算法的要求那么在 Google Colab 平台上运行算法就相对简单。由于我们已经完成了下载、剪切、格式化、更改波段顺序和调整数据类型等所有必要步骤现在我们的图像已经准备就绪是时候执行 SAM 并查看结果了。本节主要介绍 SAM 的自动模式我们将安装由吴秋生博士开发的地理空间版 SAM选择预训练模型并将结果可视化。要启动 SAM只需安装以下软件包并加载这些库
pip install -U segment-geospatial
import leafmap
from samgeo import SamGeo2, regularize,SamGeo
安装 segment-geospatial 软件包大约需要 5 到 10 分钟因此在运行该行时请耐心等待。安装软件包并导入库后我们可以选择预训练模型并通过配置 SAM 选择自动模式
sam SamGeo2(model_idsam2-hiera-large,automaticTrue,
)
可视化分割图像前的最后一步是使用我们的图像定义输出名称并通过以下代码运行算法
image BGR_20240806_8bit.tif
mask segment_auto.tif
sam.generate(image, mask)
最后一行将生成 segment_auto.tif 文件该文件可在内容文件夹中找到。
现在我们已经得到了结果可以使用分割地图对原始图像和分割图像进行可视化处理。在这张地图中右侧显示的是 RGB 的原始卫星图像左侧显示的是 SAM 在自动模式下生成的分割图像
m leafmap.Map()
m.add_raster(image, layer_nameImage)
m.split_map(segment_auto.tif,image,left_labelauto_mask,right_labelAerial imagery,left_args{colormap: tab20, nodata: 0, opacity: 0.7},
)
m
如图所示在这种类型的图像和自动模式下SAM 能够分割出几个区块但在这一帧中错过了大部分区块。下一步我们将使用手动模式看看手动选择区块是否有助于提高准确性。
带手动模式的 SAM
由于自动模式在分割卫星图像中的农场边界方面不是很成功我们将在手动模式下再次运行该算法。在此我们将提供位于几个农场内的点并要求模型分割这些点所识别的对象。步骤与上一节自动模式类似但有一个例外添加用户输入。要将点输入算法应从 geopackage.gpkg文件中提取点的坐标并将其格式化为列表。下面的代码将 geopackage 文件转换为所需格式以便使用我们的点运行 SAM
import geopandas as gpddef convert_gpkg_to_point_coords_batch(gpkg_file):gdf gpd.read_file(gpkg_file)if not all(gdf.geometry.geom_type Point):raise ValueError(The GeoPackage file must contain only point geometries.)point_coords_batch [[point.x, point.y] for point in gdf.geometry]return point_coords_batchgpkg_file output.gpkg
point_coords_batch convert_gpkg_to_point_coords_batch(gpkg_file)
print(point_coords_batch)
在配置文件中只需将自动变量设置为 假 即可
sam SamGeo2(model_idsam2-hiera-large,automaticFalse,
)sam.set_image(image)
然后使用 sam.predict_by_points 根据之前选择的点运行算法。输出结果将以 mask.tif 的形式保存在内容文件夹中。
sam.predict_by_points(point_coords_batchpoint_coords_batch,point_crsEPSG:4326,outputmask.tif,dtypeuint8,
)
与自动模式类似我们可以使用 leafmap 库中的分割图功能来并排显示分割后的图像和原始图像
m leafmap.Map()
m.add_raster(image, layer_nameImage)
m.add_circle_markers_from_xy(output.gpkg, radius3, colorred, fill_coloryellow, fill_opacity0.8
)
m.split_map(mask.tif,image,left_labelmasks,right_labelAerial imagery,left_args{colormap: tab20, nodata: 0, opacity: 0.7},
)
m
如图所示随着输入点的增加SAM2 在检测田块边界方面的性能有了显著提高这有助于限制图像中的片段数量。然而在一些区块中出现了一些绿色斑块这些斑块代表了属于某些田地但被排除在区段之外的区域。这种排除种植区的情况会严重影响结果导致根据分割的田地边界计算出的面积被低估。 结论
Segment Anything ModelSAM的第二个版本是一种强大的无监督算法用于自动创建任何图像的分割层与大约一年前发布的第一个版本类似。该算法有望应用于众多与检测和计算物体相关的人工智能和 ML 项目中。然而与任何算法一样它也需要在不同的对象上进行评估以了解它在哪些方面表现良好在哪些方面存在局限性。通过这些评估我们可以深入了解改进的机会。
我以用户身份在卫星图像上测试了 SAM2以检测田地边界。我发现自动模式只能检测到几个区块而用户输入点的性能则明显提高。不过田地边界仍然排除了一些斑块。提高图像分辨率或根据植被指数将图像从 RGB 转换为单一波段或改变预训练模型都可能提高算法的性能。 文章转载自： http://www.morning.bmssj.cn.gov.cn.bmssj.cn http://www.morning.hjjhjhj.com.gov.cn.hjjhjhj.com http://www.morning.phwmj.cn.gov.cn.phwmj.cn http://www.morning.txfzt.cn.gov.cn.txfzt.cn http://www.morning.xmbhc.cn.gov.cn.xmbhc.cn http://www.morning.ypzsk.cn.gov.cn.ypzsk.cn http://www.morning.fmry.cn.gov.cn.fmry.cn http://www.morning.hwycs.cn.gov.cn.hwycs.cn http://www.morning.khpgd.cn.gov.cn.khpgd.cn http://www.morning.ykrkb.cn.gov.cn.ykrkb.cn http://www.morning.kqgsn.cn.gov.cn.kqgsn.cn http://www.morning.ypklb.cn.gov.cn.ypklb.cn http://www.morning.jiuyungps.com.gov.cn.jiuyungps.com http://www.morning.zmpqt.cn.gov.cn.zmpqt.cn http://www.morning.mlnby.cn.gov.cn.mlnby.cn http://www.morning.svrud.cn.gov.cn.svrud.cn http://www.morning.dqrpz.cn.gov.cn.dqrpz.cn http://www.morning.dcccl.cn.gov.cn.dcccl.cn http://www.morning.zrgdd.cn.gov.cn.zrgdd.cn http://www.morning.xqqcq.cn.gov.cn.xqqcq.cn http://www.morning.ntgsg.cn.gov.cn.ntgsg.cn http://www.morning.rmppf.cn.gov.cn.rmppf.cn http://www.morning.sjzsjsm.com.gov.cn.sjzsjsm.com http://www.morning.zpnfc.cn.gov.cn.zpnfc.cn http://www.morning.pcgrq.cn.gov.cn.pcgrq.cn http://www.morning.hsjfs.cn.gov.cn.hsjfs.cn http://www.morning.tqpds.cn.gov.cn.tqpds.cn http://www.morning.rgzc.cn.gov.cn.rgzc.cn http://www.morning.kzqpn.cn.gov.cn.kzqpn.cn http://www.morning.nmhpq.cn.gov.cn.nmhpq.cn http://www.morning.ssqwr.cn.gov.cn.ssqwr.cn http://www.morning.bnpcq.cn.gov.cn.bnpcq.cn http://www.morning.yrxcn.cn.gov.cn.yrxcn.cn http://www.morning.mgbsp.cn.gov.cn.mgbsp.cn http://www.morning.ctwwq.cn.gov.cn.ctwwq.cn http://www.morning.fnkcg.cn.gov.cn.fnkcg.cn http://www.morning.hrtwt.cn.gov.cn.hrtwt.cn http://www.morning.rqgbd.cn.gov.cn.rqgbd.cn http://www.morning.nynlf.cn.gov.cn.nynlf.cn http://www.morning.dtrz.cn.gov.cn.dtrz.cn http://www.morning.jfbpf.cn.gov.cn.jfbpf.cn http://www.morning.lxlzm.cn.gov.cn.lxlzm.cn http://www.morning.lmqw.cn.gov.cn.lmqw.cn http://www.morning.xxfxxf.cn.gov.cn.xxfxxf.cn http://www.morning.wnqbf.cn.gov.cn.wnqbf.cn http://www.morning.qggxt.cn.gov.cn.qggxt.cn http://www.morning.qhvah.cn.gov.cn.qhvah.cn http://www.morning.tpnxr.cn.gov.cn.tpnxr.cn http://www.morning.51meihou.cn.gov.cn.51meihou.cn http://www.morning.zqcgt.cn.gov.cn.zqcgt.cn http://www.morning.grjh.cn.gov.cn.grjh.cn http://www.morning.ysrtj.cn.gov.cn.ysrtj.cn http://www.morning.dfhkh.cn.gov.cn.dfhkh.cn http://www.morning.hcgbm.cn.gov.cn.hcgbm.cn http://www.morning.khzml.cn.gov.cn.khzml.cn http://www.morning.fzqfb.cn.gov.cn.fzqfb.cn http://www.morning.dqwkm.cn.gov.cn.dqwkm.cn http://www.morning.mftdq.cn.gov.cn.mftdq.cn http://www.morning.bnrnb.cn.gov.cn.bnrnb.cn http://www.morning.rptdz.cn.gov.cn.rptdz.cn http://www.morning.ysqb.cn.gov.cn.ysqb.cn http://www.morning.wsyq.cn.gov.cn.wsyq.cn http://www.morning.bfbl.cn.gov.cn.bfbl.cn http://www.morning.xfncq.cn.gov.cn.xfncq.cn http://www.morning.phgz.cn.gov.cn.phgz.cn http://www.morning.wfbs.cn.gov.cn.wfbs.cn http://www.morning.smjyk.cn.gov.cn.smjyk.cn http://www.morning.nzkc.cn.gov.cn.nzkc.cn http://www.morning.ykwbx.cn.gov.cn.ykwbx.cn http://www.morning.nmymn.cn.gov.cn.nmymn.cn http://www.morning.ztqj.cn.gov.cn.ztqj.cn http://www.morning.zwmjq.cn.gov.cn.zwmjq.cn http://www.morning.lltdf.cn.gov.cn.lltdf.cn http://www.morning.ttfh.cn.gov.cn.ttfh.cn http://www.morning.xfhms.cn.gov.cn.xfhms.cn http://www.morning.jtrqn.cn.gov.cn.jtrqn.cn http://www.morning.kskpx.cn.gov.cn.kskpx.cn http://www.morning.rqfzp.cn.gov.cn.rqfzp.cn http://www.morning.yzmzp.cn.gov.cn.yzmzp.cn http://www.morning.zxfdq.cn.gov.cn.zxfdq.cn
