遥感影像计算Jaccard 相似性系数

遥感影像可以为矢量和栅格,这里以栅格为例。影响格式为二值图像

代码语言:javascript
复制
#import rasterio and sklearn
import rasterio
from rasterio.warp import reproject, Resampling, calculate_default_transform
from sklearn.metrics import jaccard_score

#read the first raster image
with rasterio.open("G:\aa\二值图像\globeland301.tif") as src1:
raster1 = src1.read(1)
profile1 = src1.profile

#read the second raster image
with rasterio.open("G:\aa\二值图像\CGLS-LC1001.tif") as src2:
raster2 = src2.read(1)
profile2 = src2.profile

#resample the first raster to 100m resolution
resampled_raster1 = rasterio.open("resampled_raster1.tif", "w", **profile1)
resampled_raster1.write(raster1, 1)
resampled_raster1.resample(100, Resampling.nearest)
resampled_raster1.close()

#resample the second raster to 100m resolution
resampled_raster2 = rasterio.open("resampled_raster2.tif", "w", **profile2)
resampled_raster2.write(raster2, 1)
resampled_raster2.resample(100, Resampling.nearest)
resampled_raster2.close()

#calculate the destination transform and crs for reprojecting
dst_transform, dst_width, dst_height = calculate_default_transform(
resampled_raster1.crs, "EPSG:4326",
resampled_raster1.width, resampled_raster1.height,
*resampled_raster1.bounds)

#reproject the first raster to WGS-84
reprojected_raster1 = rasterio.open("reprojected_raster1.tif", "w", crs="EPSG:4326",
transform=dst_transform,
width=dst_width,
height=dst_height,
nodata=0) #set nodata value to 0
reproject(
source=rasterio.band(resampled_raster1, 1),
destination=rasterio.band(reprojected_raster1, 1),
src_transform=resampled_raster1.transform,
src_crs=resampled_raster1.crs,
dst_transform=reprojected_raster1.transform,
dst_crs=reprojected_raster1.crs,
resampling=Resampling.nearest)
reprojected_raster1.close()

#reproject the second raster to WGS-84
reprojected_raster2 = rasterio.open("reprojected_raster2.tif", "w", crs="EPSG:4326",
transform=dst_transform,
width=dst_width,
height=dst_height,
nodata=0) #set nodata value to 0
reproject(
source=rasterio.band(resampled_raster2, 1),
destination=rasterio.band(reprojected_raster2, 1),
src_transform=resampled_raster2.transform,
src_crs=resampled_raster2.crs,
dst_transform=reprojected_raster2.transform,
dst_crs=reprojected_raster2.crs,
resampling=Resampling.nearest)
reprojected_raster2.close()

#calculate Jaccard similarity coefficient between the two rasters
j_raster = jaccard_score(reprojected_raster1.read(1).flatten(), reprojected_raster2.read(1).flatten())

#print the result
print("Jaccard similarity coefficient between the two rasters:", j_raster)