alex_bn_lee

导航

【854】通过polygon切取tif栅格数据(rasterio读取与存储)

参考:Cutting a polygon from TIFF with Python [closed]


import rasterio
import rasterio.mask
import geopandas as gpd

dataset = rasterio.open("wc2.1_10m_elev.tif")
gdf_africa = gpd.read_file("africa1_map.gpkg")

poly = gdf_africa.loc[0, "geometry"]

# out_image is the tiff image cut by the poly
out_image, out_transform = rasterio.mask.mask(dataset, [poly], crop=True)

下面的例子是两幅世界地图,一幅是矢量一幅是栅格,通过矢量把栅格剪切,并获取其平均值!

polys = list(gdf_africa["geometry"]) 
avg_values = []

for poly in polys:
    out_image, out_transform = rasterio.mask.mask(dataset, [poly], crop=True)

    out_image[out_image == -32768] = 0
    sum_values = out_image.sum()
    out_image[out_image != 0] = 1 
    num_values = out_image.sum()

    avg_values.append(sum_values/num_values) 

-32768对应的为最小值,也就是非数值的部分,需要去掉,之后取和,就是总体象元值的和。

再把非0的数据赋值为1,取和之后就是计算有值的象元的个数,然后两者取商就是平均值。

这样计算效率很高,最开始我是考虑遍历整个地图,然后通过判断像素点是否在polygon里面计算,效率非常低!

存储为 tiff 文件

import rasterio
import rasterio.mask
import geopandas as gpd


# Read roofs
gdf = gpd.read_file("input_shapes.geojson")  # Your roofs
gdf.set_index("roof_id")
roof = gdf.loc[1234]  # Your roof id


# Open input raster and write masked (clipped) output raster
with rasterio.open("input_raster.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, [roof["geometry"]], crop=True)
    out_meta = src.meta

    out_meta.update(
        {
            "height": out_image.shape[1],
            "width": out_image.shape[2],
            "transform": out_transform,
        }
    )

    with rasterio.open("output_raster.tif", "w", **out_meta) as dest:
        dest.write(out_image)

 

posted on 2023-07-06 08:25  McDelfino  阅读(409)  评论(0)    收藏  举报