基于python制作陆地掩膜

想要一个和输入raster大小和投影都一样的掩膜，然后找到了以下代码，这里的110就是给陆地付的值。速度算很快

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 23 16:41:35 2024

@author: Asus
"""
import numpy as np
from osgeo import gdal
import rasterio
from rasterio.mask import mask
from geopandas import GeoDataFrame
from rasterio import features

def extract_by_shp_rasterio(in_shp_path, in_raster_path, out_raster_path):
    maskdata = GeoDataFrame.from_file(in_shp_path)
    # with fiona.open(in_shp_path, "r", encoding='utf-8') as shapefile:
    #     # 获取所有要素feature的形状geometry
    #     geoms = [feature["geometry"] for feature in shapefile]
    
    
    # 裁剪陆地的,即去除shp里的
    with rasterio.open(in_raster_path) as src:
        shpdata = maskdata.to_crs(src.crs)  # 投影转换，使矢量数据与栅格数据投影参数一致
        geoms = [feature for feature in shpdata.geometry]
        out_image, out_transform = mask(src, geoms, invert=True, nodata=110,all_touched=True)
        out_meta = src.meta.copy()
        # before = src.read(1)
        # out_image[before[None]==3]=3
    
    
    
    
    # 更新输出文件的元数据
    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],
                     "width": out_image.shape[2],
                     "transform": out_transform})
    
    # 保存
    with rasterio.open(out_raster_path, "w", **out_meta) as dest:
        dest.write(out_image)
        

    
if __name__ == "__main__":
    in_shp_path=r"G:\mask\simple\simple_land_1kmBuffer.shp"
    inputpath=r"Z:\download\merge\monthly\frequency\201601_frequency.tif"
    outputpath=r"G:\simple\simple_land_1kmBuffer_Arctic40mlandmask110.tif"
    extract_by_shp_rasterio(in_shp_path, inputpath, outputpath)