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

淡淡的烟草味﹌ 2023-10-13 22:45 14阅读 0赞

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

#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)