想要一个和输入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)