def clip_veg(cities, raster_path, row_index):
with rasterio.open(raster_path) as src:
city = cities.iloc[row_index]
city_geometry = city['geometry']
out_image, out_transform = mask(src, [city_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
})
print(out_image.shape)
return out_image, out_transform, out_meta
def veg_topoints(image, transform):
rows, cols = np.where((image > 0) & (image < 255))
x_coords = transform[0] + (cols + 0.5) * transform[1]
y_coords = transform[3] + (rows + 0.5) * transform[5]
geometries = [Point(x, y) for x, y in zip(x_coords, y_coords)]
values = image[rows, cols]
return gpd.GeoDataFrame({'raster_val': values, 'geometry': geometries}, crs='epsg:4326')
写了两个函数对数据依次操作,之后结果竟然报错,image是三维数组,查了一下rasterio.mask.mask的用法,裁剪之后是三维数组,就算原数据只有一个波段,最后也是(1,height, width),所以下一个函数np.where()部分就报错了。加入image=image[0],变成二维数组就好了。