之前搜索了rasterio的矢量裁剪的方法,但是真正运用的时候报错ValueError: No valid geometry objects found for rasterize,但是实际上我的shp和raster是重叠的,发现是没有将矢量数据的投影与栅格的进行统一。因此采用下面的方法即可:
import fiona
import rasterio
from rasterio.mask import mask
from geopandas import GeoDataFrame
from rasterio import features
from shapely.geometry import shape
from scipy import ndimage
try:
import gdal
except:
from osgeo import gdal
def extract_by_shp_rasterio(in_shp_path, in_raster_path, out_raster_path):
maskdata = GeoDataFrame.from_file(in_shp_path)
# crop即保留shp里的,invert=True即去除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, crop=True, nodata=4,all_touched=True)
out_meta = src.meta.copy()
# 更新输出文件的元数据
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)
另外的一个需求是需要筛选shp的属性,比如说本次我需要FP这个属性=='80',因此代码做了改动,注意直接shpdata.属性名即可获取属性
def extract_by_landfastshp_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]
# 裁剪
with rasterio.open(in_raster_path) as src:
shpdata = maskdata.to_crs(src.crs) # 投影转换,使矢量数据与栅格数据投影参数一致
# geoms = [feature for feature in shpdata if feature.FP=='80']
geoms=[]
for i in range(0, len(shpdata) ):
geo = shpdata.geometry[i] #获取空间属性,即GeoSeries
name = shpdata.FP[i] #获得属性数据,即字段
if name =='08':
geoms.append(geo)
# geoms = [feature for feature in shpdata.geometry if feature.FP=='80']
out_image, out_transform = mask(src, geoms, invert=True, nodata=6,all_touched=True)
out_meta = src.meta.copy()
# 更新输出文件的元数据
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)