【Python&RS】Pyrsgis库安装+基础函数使用教程

​

        pyrsgis库是一个用于处理地理信息系统(GIS)数据的Python库。它提供了一组功能强大的工具，可以帮助开发人员使用Python语言创建、处理、分析和可视化GIS数据。通过使用pyrsgis库，开发人员可以更轻松地理解和利用地理信息。

        pyrsgis库包含了许多常见的GIS操作和功能，例如读取和写入shapefile文件、转换坐标系、执行空间查询、计算地理特征属性等。它提供了许多方便使用的类和方法，例如GeoPandas、Shapely、Fiona、Rasterio、Pyproj和GDAL等，这些都可以帮助开发人员更高效地处理GIS数据。

一、Pyrsgis库安装

        Pyrsgis可以直接通过pip install pyrsgis安装，同样也可以下载压缩包然后本地安装。PyPI中Pyrsgis包下载地址：pyrsgis · PyPI

​

二、导入库和函数

        这些都是我后面代码需要使用到的函数，注意要导入，别到时候报错。

import os
from pyrsgis import raster, convert, ml

三、基础操作代码展示

1）获取影像基本信息

def Get_data(filepath):
    # 获取影像基本信息
    ds, data_arr = raster.read(filepath)  # 基础信息资源和数组
    ds_bands = ds.RasterCount  # 波段数
    ds_width = ds.RasterXSize  # 宽度
    ds_height = ds.RasterYSize  # 高度
    ds_bounds = ds.bbox  # 四至范围
    ds_geo = ds.GeoTransform  # 仿射地理变换参数
    ds_prj = ds.Projection  # 投影坐标系
    print("影像的宽度为：" + str(ds_width))
    print("影像的高度为：" + str(ds_height))
    print("仿射地理变换参数为：" + str(ds_geo))
    print("投影坐标系为：" + str(ds_prj))

2）计算NDVI

        这里给大家介绍一个经典案例，就是NDVI的计算。通过这个应该很容易就能理解Pyrsgis库的数据结构了。

def Get_NDVI(filepath):
    # 计算NDVI
    ds, data_arr = raster.read(filepath)
    red_arr = data_arr[3, :, :]
    nir_arr = data_arr[4, :, :]
    result_arr = (nir_arr - red_arr) / (nir_arr + red_arr)
    # result_arr = (data_arr[4, :, :] - data_arr[3, :, :]) / (data_arr[4, :, :] + data_arr[3, :, :])
    output_file = r'E:/path_to_your_file/landsat8_result.tif'
    raster.export(result_arr, ds, output_file, dtype='float32', bands="all", nodata=0, compress="LZW")
    # 写入的数组，基础信息，路径，格式，波段，无效值，压缩方式

3）空间位置裁剪

        这里的裁剪主要是按照输入的空间矩形进行裁剪，并没有演示如何使用shp进行裁剪。这个可以应用于分幅裁剪、滑动裁剪等。空行分割的是实现这个功能的两种函数的使用方式。

def Clip_data(filepath):
    # 按掩膜提取
    ds, data_arr = raster.read(filepath)
    print('Original bounding box:', ds.bbox)
    print('Original shape of raster:', data_arr.shape)
    new_ds, clipped_arr = raster.clip(ds, data_arr, x_min=770000, x_max=790000, y_min=1420000, y_max=1440000)
    raster.export(clipped_arr, new_ds, r'E:/path_to_your_file/clipped_file.tif')

    infile = r'E:/path_to_your_file/your_file.tif'
    outfile = r'E:/path_to_your_file/clipped_file.tif'
    raster.clip_file(infile, x_min=770000, x_max=790000, y_min=1420000, y_max=1440000, outfile=outfile)

4）移除无效值

        这里的函数是移除无效值的，如-9999之类的，理论上应该也可以修改其他的DN值，但我自己没去试过，大家可以自行尝试。

def Modify_data(filepath):
    # 修改影像数组，如移除无效值
    ds, data_arr = raster.read(filepath)
    new_ds, new_arr = raster.trim(ds, data_arr, remove=-9999)
    print('Shape of the input array:', data_arr.shape)
    print('Shape of the trimmed array:', new_arr.shape)

    ds, data_arr = raster.read(filepath)
    new_arr = raster.trim_array(data_arr, remove=-9999)
    print('Shape of the input array:', data_arr.shape)
    print('Shape of the trimmed array:', new_arr.shape)

    infile = r'E:/path_to_your_file/your_file.tif'
    outfile = r'E:/path_to_your_file/trimmed_file.tif'
    raster.trim_file(infile, -9999, outfile)

5）平移影像

        按照x、y方向进行影像平移，可选像素和坐标进行平移。

def Shift_data(filepath):
    # 平移影像
    ds, data_arr = raster.read(filepath)
    new_ds = raster.shift(ds, x=10, y=10)  # x,y方向偏移量。按栅格的投影单位移动数据源 或分别按细胞数
    print('Original bounding box:', ds.bbox)
    print('Modified bounding box:', new_ds.bbox)

    new_ds = raster.shift(ds, x=10, y=10, shift_type='cell')  # shift_type='coordinate'
    print('Modified bounding box:', new_ds.GeoTransform)
    raster.export(data_arr, new_ds, r'E:/path_to_your_file/shifted_file.tif')

    infile = r'E:/path_to_your_file/your_file.tif'
    outfile = r'E:/path_to_your_file/shifted_file.tif'
    raster.shift_file(infile, x=10, y=10, outfile=outfile, shift_type='cell')

6）数组、表、CSV互转（包含剔除值）

        这里的函数是数组、表、CSV互转，在转换的同时可以通过参数移除某些DN值。

def Convert_data(filepath):
    # 数组转表、CSV，修改值
    input_file = r'E:/path_to_your_file/raster_file.tif'
    ds, data_arr = raster.read(input_file)  # Shape of the input array: (6, 800, 400)
    data_table = convert.array_to_table(data_arr)  # Shape of the reshaped array: (320000, 6)
    # 该函数将单波段或多波段栅格数组转换为表，其中 列表示输入波段，每行表示一个单元格。

    input_file = r'E:/path_to_your_file/raster_file.tif'
    ds, data_arr = raster.read(input_file)
    data_table = convert.array_to_table(data_arr)
    print('Shape of the input array:', data_arr.shape)  # Shape of the input array: (6, 800, 400)
    print('Shape of the reshaped array:', data_table.shape)  # Shape of the reshaped array: (320000, 6)
    input_file = r'E:/path_to_your_file/raster_file.tif'
    new_data_arr = convert.table_to_array(data_table, n_rows=ds.RasterYSize, n_cols=ds.RasterXSize)
    print('Shape of the array with newly added bands:', new_data_arr.shape)
    # Shape of the array with newly added bands: (8, 800, 400)
    new_data_arr = convert.table_to_array(data_table[:, -2:], n_rows=ds.RasterYSize, n_cols=ds.RasterXSize)
    print('Shape of the array with newly added bands:', new_data_arr.shape)
    # Shape of the array with newly added bands: (2, 800, 400)
    # 表转数组

    input_file = r'E:/path_to_your_file/raster_file.tif'
    output_file = r'E:/path_to_your_file/tabular_file.csv'
    convert.raster_to_csv(input_file, filename=output_file)
    input_dir = r'E:/path_to_your_file/'
    output_file = r'E:/path_to_your_file/tabular_file.csv'
    convert.raster_to_csv(input_dir, filename=output_file)
    convert.raster_to_csv(input_dir, filename=output_file, negative=False, remove=[10, 54, 127], badrows=False)
    # 数组转表，可剔除负值、目标值、坏波段

    input_file = r'E:/path_to_your_file/raster_file.tif'
    out_csvfile = input_file.replace('.tif', '.csv')
    convert.raster_to_csv(input_file, filename=out_csvfile, negative=False)
    new_csvfile = r'E:/path_to_your_file/predicted_file.tif'
    out_tiffile = new_csvfile.replace('.csv', '.tif')
    convert.csv_to_raster(new_csvfile, ref_raster=input_file, filename=out_tiffile, compress='DEFLATE')
    convert.csv_to_raster(new_csvfile, ref_raster=input_file, filename=out_tiffile,
                          cols=['Blue', 'Green', 'KMeans', 'RF_Class'], compress='DEFLATE')
    # 数组将堆叠并导出为多光谱文件
    convert.csv_to_raster(new_csvfile, ref_raster=input_file, filename=out_tiffile,
                          cols=['Blue', 'Green', 'KMeans', 'RF_Class'], stacked=False, compress='DEFLATE')
    # 将每列导出为单独的波段，请将参数设置为 。stacked=False

7）制作深度学习标签

        此函数根据单波段或多波段栅格阵列生成影像片。图像芯片可以用作深度学习模型的直接输入（例如。卷积神经网络），输出格式：(4198376, 7, 7, 6)

def Create_CNN(filepath):
    # 此函数根据单波段或多波段栅格阵列生成影像片。图像 芯片可以用作深度学习模型的直接输入（例如。卷积神经网络）
    # -----------------------------数组生成深度学习芯片-----------------------------
    infile = r'E:/path_to_your_file/your_file.tif'
    ds, data_arr = raster.read(infile)
    image_chips = ml.array_to_chips(data_arr, y_size=7, x_size=7)
    print('Shape of input array:', data_arr.shape)  # Shape of input array: (6, 2054, 2044)
    print('Shape of generated image chips:', image_chips.shape)  # Shape of generated image chips: (4198376, 7, 7, 6)

    infile = r'E:/path_to_your_file/your_file.tif'
    ds, data_arr = raster.read(infile)
    image_chips = ml.array2d_to_chips(data_arr, y_size=5, x_size=5)
    print('Shape of input array:', data_arr.shape)  # Shape of input array: (2054, 2044)
    print('Shape of generated image chips:', image_chips.shape)  # Shape of generated image chips: (4198376, 5, 5)

    # ----------------------------影像直接生成深度学习芯片----------------------------
    infile_2d = r'E:/path_to_your_file/your_2d_file.tif'
    image_chips = ml.raster_to_chips(infile_2d, y_size=7, x_size=7)
    print('Shape of single band generated image chips:', image_chips.shape)
    # Shape of single bandgenerated image chips: (4198376, 7, 7)

    infile_3d = r'E:/path_to_your_file/your_3d_file.tif'
    image_chips = ml.raster_to_chips(infile_3d, y_size=7, x_size=7)
    print('Shape of multiband generated image chips:', image_chips.shape)
    # Shape of multiband generated image chips: (4198376, 7, 7, 6)

8）翻转影像

        按照东西或南北方向翻转影像

def Reverse_Image(filepath):
    # 按照东西、南北方向反转影像
    # -------------------------------北向、东向翻转--------------------------------
    input_file = r'E:/path_to_your_file/your_file.tif'
    ds, data_arr = raster.read(input_file)
    north_arr, east_arr = raster.north_east(data_arr)
    print(north_arr.shape, east_arr.shape)

    north_arr, east_arr = raster.north_east(data_arr, flip_north=True, flip_east=True)
    north_arr = raster.north_east(data_arr, layer='north')
    from matplotlib import pyplot as plt
    plt.imshow(north_arr)
    plt.show()
    plt.imshow(east_arr)
    plt.show()

    input_file = r'E:/path_to_your_file/your_file.tif'
    ds, data_arr = raster.read(input_file)
    north_arr, east_arr = raster.north_east(data_arr)
    print(north_arr.shape, east_arr.shape)
    north_arr = raster.north_east(data_arr, layer='north')
    from matplotlib import pyplot as plt
    plt.imshow(north_arr)
    plt.show()
    plt.imshow(east_arr)
    plt.show()
    raster.export(north_arr, ds, r'E:/path_to_your_file/northing.tif', dtype='float32')
    raster.export(east_arr, ds, r'E:/path_to_your_file/easting.tif', dtype='float32')
    # -------------------------使用参考.tif文件生成北向栅格----------------------------
    reference_file = r'E:/path_to_your_file/your_file.tif'
    raster.northing(file1, r'E:/path_to_your_file/northing_number.tif', flip=False, value='number')
    raster.northing(file1, r'E:/path_to_your_file/northing_normalised.tif', value='normalised')  # 输出栅格进行归一化
    raster.northing(file1, r'E:/path_to_your_file/northing_coordinates.tif', value='coordinates')
    raster.northing(file1, r'E:/path_to_your_file/northing_number_compressed.tif', compress='DEFLATE')

    reference_file = r'E:/path_to_your_file/your_file.tif'
    raster.easting(file1, r'E:/path_to_your_file/easting_number.tif', flip=False, value='number')
    raster.easting(file1, r'E:/path_to_your_file/easting_normalised.tif', value='normalised')
    raster.easting(file1, r'E:/path_to_your_file/easting_normalised.tif', value='normalised')
    raster.easting(file1, r'E:/path_to_your_file/easting_number_compressed.tif', compress='DEFLATE')

四、总结

        Pyrsgis库之前使用的时候是因为要进行卷积神经网络的深度学习，然后里面制作深度学习标签的函数还是不错的，可以用一行代码实现标签的制作。但是如果数据过大，内存就会溢出报错，这个是Pyrsgis库没有解决的，当然我也没解决=。=大家可以自己尝试一下，有解决办法可以和我分享一下。总的来说Pyrsgis和Rasterio这两个库都还不错，都在GDAL的基础上进行了二开，方便了很多操作。

​