python - 无法使用 rasterio.mask 剪辑图像
问题描述
我正在尝试使用 python 中的形状或 geojson 文件来剪辑我的 tiff 文件。剪切图像的代码是 -
from datetime import date
import geopandas as gpd
import rasterio
import rasterio.features
import rasterio.warp
from shapely.geometry import MultiPolygon, Polygon
import subprocess
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.mask import mask
nReserve = gpd.read_file('poly.shp')
# the polygon GeoJSON geometry
nReserve_proj = nReserve.to_crs({'init': 'epsg:32643'})
with rasterio.open("RGBNew.tiff") as src:
out_image, out_transform = rasterio.mask.mask(src, nReserve_proj.geometry,crop=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("RGB_masked.tif", "w", **out_meta) as dest:
dest.write(out_image)
创建上述使用的形状文件的代码是 -
import pandas as pd
from shapely.geometry import Polygon
def bbox(lat,lng, margin):
#return (Polygon([[19.386508042365318,72.83352971076965],[19.386163944316234,72.83352971076965],[19.387256959134895,72.83352971076965],[19.38746948894201,72.83352971076965],[19.386508042365318,72.83352971076965]]))
return (Polygon([[72.83352971076965,19.386508042365318],[72.83352971076965,19.386163944316234],[72.83352971076965,19.387256959134895],[72.83352971076965,19.38746948894201],[72.83352971076965,19.386508042365318]]))
gpd.GeoDataFrame(crs = {'init':'epsg:32643'},geometry = [bbox(10,10, 0.25)]).to_file('poly.shp')
但我得到了错误-
文件“tryWithNewEPSG.py”,第 24 行,在 out_image,out_transform = mask(src, geoms,crop=True) 文件“/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py ",第 181 行,在掩码 pad=pad) 文件 "/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py",第 87 行,在 raster_geometry_mask raise ValueError('输入形状做不重叠栅格。')ValueError:输入形状不重叠栅格。
我目前正在使用代码检查 epsg-
import rasterio
with rasterio.open('NDVI.tif') as src:
print (src.crs)
并已确认相同;我什至尝试将两者的 epsg 更改为 4326;但仍然没有工作。
out_image, out_transform = rasterio.mask.mask(src, nReserve_proj.geometry,crop=True)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
解决方案
问题是创建形状文件的方法;使用以下代码创建形状文件-
def create_polygon(coords):
ring = ogr.Geometry(ogr.wkbLinearRing)
for coord in coords:
ring.AddPoint(coord[0], coord[1])
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
return poly.ExportToWkt()
def write_shapefile(poly, out_shp):
"""
https://gis.stackexchange.com/a/52708/8104
"""
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(out_shp)
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
## If there are multiple geometries, put the "for" loop here
# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', 123)
# Make a geometry, from Shapely object
#geom = ogr.CreateGeometryFromWkb(poly.wkb)
geom = ogr.CreateGeometryFromWkt(poly)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
# Save and close everything
ds = layer = feat = geom = None
def main(coords, out_shp):
poly = create_polygon(coords)
write_shapefile(poly, out_shp)
if __name__ == "__main__":
#coords = [(73.0367660522461,18.979999927217715), (73.03436279296875,18.96019467106038), (73.05976867675781,18.96214283338193), (73.06011199,18.979999927217715),(73.0367660522461,18.979999927217715)]
coords = [(73.97164740781432,18.527929607251075),(73.97185125569945,18.528550143773767),(73.97234478215819,18.528305998525394),(73.97215166310912,18.52775667044178),(73.97164740781432,18.527929607251075)]
out_shp = r'polygon.shp'
main(coords, out_shp)
在坐标中提到坐标;确保您正在剪辑的 tiff 具有epsg:4326您可以使用转换代码来执行此操作-
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
dst_crs = 'epsg:4326'
with rasterio.open('NDVIData.tiff') as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open('NDVINew.tiff', 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
其中 NDVIData.tiff 是我的原始 tiff,NDVINew.tiff 是我的新 tiff,带有新的 epsg 4326。
现在尝试运行剪辑代码
推荐阅读
- android - 安卓应用证书
- python - “TypeError:如果启用引用,则必须设置quotechar”在csv模块中
- flutter - Flutter 转换货币格式
- c# - 参考字符串参数为空?
- string - 如何将字符串列转换为具有浮点值pyspark的udt向量
- css - 移动视图上的重叠文本
- javascript - React 在按钮单击时渲染新图像(每次检索图像的新 fetch 调用)
- amazon-web-services - 使用 AWS 负载均衡器之一对 gRPC 请求进行负载均衡
- swift - 如何禁用允许删除提示?
- reactjs - 静态和动态Key和Value的setState的区别?