Geopandas + rasterio:将矢量隔离为 png

Geopandas + rasterio : isolate a vector as png

我正在尝试将城市边界隔离为 png 的单个部分。我的目标是将此 png 叠加到非常旧的卫星照片上。

为此,我收集了一个复制照片尺寸的光栅文件和一个带边界的矢量文件。然后,我使用光栅:

import rasterio
from rasterio.plot import show

src = rasterio.open("my_raster.tiff")

与 geopandas 的等价物:

import geopandas as gpd

GDF = gpd.read_file("boundary.shp")

我检查了 src 和 GDF 之间的坐标参考系统完全相同,然后我使用 matplotlib 正确放置了边界:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20, 10))
show(src.read(), transform=src.transform, ax=ax)
GDF.plot(ax=ax, color='white')
plt.show()

其中显示:

效果很好,但我无法使用 savefig() 仅保存 png 中的边界。我试图将 ax 分开,ax1 用于光栅,ax2 用于矢量,但它没有用...

我可以只保存这部分图吗?

好的,经过一些搜索,我将我的情节保存如下:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20, 10))
show(src.read(), transform=src.transform, ax=ax)
GDF.plot(ax=ax, color='#fff')
ax.set_axis_off()
fig.savefig("test.png", dpi=220, bbox_inches = 'tight')
plt.show()

然后,我在一个小函数中使用了 PIL :

from PIL import ImageTk, Image

def only_boundary(image):
    # first, convert picture as RGBA
    with Image.open(image).convert("RGBA") as img:
        pixels = img.load()
        for i in range(img.size[0]):
            for j in range(img.size[1]):
                # if a pixel is not white...
                if pixels[i,j] != (255, 255, 255,255):
                    #it becomes transparent
                    pixels[i,j] = (0, 0, 0, 0)
    # then the loops are over, we save
    im = img.save(image)

而only_boundary("test.png")保存了好结果!

我们可以使用 OpenCV 等高线使方法更清晰。

myfilter='example'
myfilter_raster=os.path.join(raster_path,myfilter+'.tif')
with rasterio.open(myfilter_raster) as src:    
vector_df=gdf_rbb[gdf_rbb.idarpt==myfilter].copy()

out_image, out_transform = rasterio.mask.mask(src, vector_df.geometry.to_list(), crop=False)
out_meta = src.meta        

out_meta.update({"driver": "PNG",
              "height": out_image.shape[1],
              "width": out_image.shape[2],
              "transform": out_transform})

out_file=myfilter+'.png'
out_file=os.path.join(mask_images_path,out_file)
print('Generated' ,out_file)

mask =  out_image[0].astype("uint8")
mask[mask > 0] = 255

border = cv2.copyMakeBorder(mask, 1, 1, 1, 1, cv2.BORDER_CONSTANT, value=0 )
contours, hierarchy = cv2.findContours(border, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE, offset=(-1, -1))

boundary_image=np.zeros(mask.shape)
for contour in contours:
    cv2.drawContours(boundary_image,[contour],0,(255,255,255),3)
plt.imshow(boundary_image)
plt.show()

with rasterio.open(out_file, 'w', **out_meta) as dst:
    dst.write(boundary_image , 1)

# from rasterio.plot import show
# fig, ax = plt.subplots(figsize=(20, 10))
# show(src.read(), transform=src.transform, ax=ax)
# vector_df.plot(ax=ax, color='white')
# plt.show()

另一种使用 opencv 边缘检测的简单方法os。

myfilter='example'
myfilter_raster=os.path.join(raster_path,myfilter+'.tif')
with rasterio.open(myfilter_raster) as src:    
    vector_df=gdf_rbb[gdf_rbb.idarpt==myfilter].copy()
    
    out_image, out_transform = rasterio.mask.mask(src, vector_df.geometry.to_list(), crop=False)
    out_meta = src.meta        
   
    out_meta.update({"driver": "PNG",
                  "height": out_image.shape[1],
                  "width": out_image.shape[2],
                  "transform": out_transform})
    
    out_file=myfilter+'2.png'
    out_file=os.path.join(mask_images_path,out_file)
    print('Generated' ,out_file)

    mask =  out_image[0].astype("uint8")
    mask[mask > 0] = 255

    edges = cv2.Canny(mask,100,200)
    plt.subplot(121)
    plt.axis('off')
    plt.imshow(mask,cmap = 'gray')
    plt.title('Original Image')
    plt.subplot(122),
    plt.imshow(edges,cmap = 'gray')
    plt.title('Edge Image')
    plt.show()
    
    with rasterio.open(out_file, 'w', **out_meta) as dst:
        dst.write(edges , 1)