rasterio 库：raster 文件操作

rasterio 库：raster 文件操作

1. 读写

import rasterio
# 
path = r'land_cover/land_cover.tif'
# 读取 zip 中的文件
path = r'zip+file:data.zip!/land_cover/land_cover.tif'
data = rasterio.open(path)

常用属性：

• .count: bound count

• .indexes: bound indices

• .width and .height: number of pixels of width and height

• .bounds: boundary

• .crs: CRS

  • 返回为 rasterio.crs.CRS 类

  • 注意：GeoDataFrame.crs 返回为 pyproj.crs.crs.CRS 类，也可以使用 equal 方法进行比较

• .res: Resolution, in CRS Unit

• .nodata: nodata value

gdf.crs.equals(raster.crs)
# out: True or False
raster.crs.equals(gdf.crs)
# AttributeError: 'rasterio.crs.CRS' object has no attribute 'equals'

常用方法：

• .read(): return numpy array

常用函数

点击查看代码

def raster_sample_point(raster, geodata):
    '''
    Parameters:
    ----
    raster (rasterio.io.DatasetReader)
    geodata (geopandas.GeoDataFrame or geopandas.GeoSeries, Point)
    
    Return (pandas.Series):
    ---
    res
    '''
    if not isinstance(geodata, gpd.GeoDataFrame):
        print('type:', type(geodata))
        return None
    
    try:
        assert geodata.crs.equals(raster.crs)
    except:
        print('raster crs is', raster.crs, 'while geodata crs is', geodata.crs)
        return None
    
    coord_list = [(x, y) for x, y in zip(geodata['geometry'].x, geodata['geometry'].y)]
    
    res = pd.Series([x[0] for x in raster.sample(coord_list)], 
                    index = geodata.index)
    return res
# ==============================================================================
def raster_buffer_sum(raster, geodata, nodata=None):
    '''
    Parameters:
    ----
    raster (rasterio.io.DatasetReader)
    geodata (geopandas.GeoDataFrame or geopandas.GeoSeries)
    
    Return:
    ---
    res (pandas.Series): using "geodata"'s index
    '''
    try:
        assert geodata.crs.equals(raster.crs)
    except:
        print('raster crs is', raster.crs, 'while geodata crs is', geodata.crs)
        return None
    
    if nodata is None:
        nodata = raster.nodata

    res = {}

    for ix, row in geodata.iterrows():
        # mask
        clip, _ = rasterio.mask.mask(raster, [row['geometry']], crop=True)
        clip[clip == nodata] = np.nan
        # print(clip)
        res[ix] = np.nansum(clip)
    
    res = pd.Series(res)

    return res
# ==============================================================================
def raster_buffer_category_count(raster, geodata, nodata=None):
    '''
    Parameters:
    ----
    raster (rasterio.io.DatasetReader)
    geodata (geopandas.GeoDataFrame or geopandas.GeoSeries)
    
    Return (pandas.DataFrame):
    ---
    res: using "geodata"'s index
    '''
    try:
        assert geodata.crs.equals(raster.crs)
    except:
        print('raster crs is', raster.crs, 'while geodata crs is', geodata.crs)
        return None
    
    from collections import Counter

    if nodata is None:
        nodata = raster.nodata
    
    res = {}

    for ix, row in geodata.iterrows():
        # mask
        clip, _ = rasterio.mask.mask(raster, [row['geometry']], crop=True)

        clip = clip.flatten().tolist()
        clip = list(filter(lambda x : x != nodata, clip))
        
        res[ix] = dict(Counter(clip))

    res = pd.DataFrame.from_dict(res, orient='index') \
            .fillna(0) \
            .astype(int) \
            .sort_index(axis='columns') 

    return res
# ==============================================================================