Search code examples
pythonmatplotlibgeopandasrasterio

Geopandas + rasterio : isolate a vector as png

I'm trying to isolate the boundary of a city as a single part of a png. My goal is to superimpose this png to very old satellite photos.

To do this, I collected a raster file which copies the dimensions of the photos and a vector file with boundary. Then, I used rasterio :

import rasterio
from rasterio.plot import show

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

And the equivalent with geopandas :

import geopandas as gpd

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

I checked that Coordinate Reference System were exactly the same between src and GDF, and then I used matplotlib to correctly put the boundary :

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

Which showed :

Rasterio + vector boundary

That worked fine, but I couldn't save only the boundary in a png with savefig(). I tried to separate ax, with an ax1 for the raster and an ax2 for the vector, but it didn't work...

Can I save only this part of the figure ?

Solution

  • We can make the approach crispier by using OpenCV contours.

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

    Yet another simple approach os to use opencv edge detection.

    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)