Python之shp文件

1、读取shp

pip install pyshp
           

pyshp(Python Shapefile Library)

是一个Python库，用于在Python脚本中对ArcGIS中的Shapefile文件(.shp，.shx，.dbf等格式)进行读写操作。

Reader类, 对shapefile文件读取；Editor类，对shapefile文件编辑；Writer类，对shapefile文件写操作

每个文件包含 “几何数据”(Geometry)和"属性数据"(Attribute Record) ,两个文件数据一一对应。

每个几何对象包含有4个属性：数据类型(shapeType)，代表该"几何数据"对象的数据类型(点，shapeType=1，线，shapeType=3，多边形，shapeType=5)；数据范围(bbox)，只针对多点数据，代表该"几何数据"对象的边界范围；数据块(parts)，只针对线或者多边形，代表该"几何数据"对象各个块的第一个点的索引；点集(points)，代表该"几何数据"对象的所有点坐标。 "属性数据"即每个"几何数据"对象在属性表中的对应项。

import shapefile
import pandas as pd

    sf = shapefile.Reader("E:/test/wanzhou_yudong_bdy.shp")
    shapes = sf.shapes()
    arr = []
    for i in range(0, len(shapes)):
        arr.append(shapes[i].bbox)
    grids = pd.DataFrame(arr, columns=['min_lon', 'min_lat', 'max_lon', 'max_lat'])
    min_lon = grids['min_lon'].min()
    min_lat = grids['min_lat'].min()
    max_lon = grids['max_lon'].max()
    max_lat = grids['max_lat'].max()
           

shp的最大最小经纬度，和上面结果一样

import shapefile

sf = shapefile.Reader("E:/tif/shp/1.shp")
min_x, min_y, max_x, max_y = sf.bbox
           

2、判断某个点是否在shp中

import shapefile
import shapely.geometry as geometry

		lons, lats = [], []
		.... # lons = np.linspace(113.77, 114.59, 50)
			 # lats = np.linspace(22.47, 22.83, 25)
        grid_lon, grid_lat = np.meshgrid(lons, lats)
        flat_lon = grid_lon.flatten()
        flat_lat = grid_lat.flatten()
        flat_points = np.column_stack((flat_lon, flat_lat))
        in_shape_points = []
        sf = shapefile.Reader("E:/test/中小河流.shp", encoding='gbk')
        shapes = sf.shapes()
        for pt in flat_points:
            # make a point and see if it's in the polygon
            if geometry.Point(pt).within(geometry.shape(shapes[0])):
                in_shape_points.append(pt)
                print("The point is in shp")
            else:
                print("The point is not in shp")
        print(in_shape_points)
           

3、gdal生成shp

import osgeo.ogr as ogr
import osgeo.osr as osr

def run():
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource("Gooise.shp")

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(32631)

    layer = data_source.CreateLayer("Gooise", srs, ogr.wkbMultiPolygon)

    # field_name = ogr.FieldDefn("Name", ogr.OFTString)
    # field_name.SetWidth(24)
    # layer.CreateField(field_name)
    #
    # field_name = ogr.FieldDefn("Car_per_house", ogr.OFTString)
    # field_name.SetWidth(24)
    # layer.CreateField(field_name)

    feature = ogr.Feature(layer.GetLayerDefn())
    # feature.SetField("Name", 'Gooise')
    # feature = ogr.Feature(layer.GetLayerDefn())
    # feature.SetField("Car_per_house", '1.2')

    wkt = 'polygon((646080 5797460,648640 5797460,648640 5794900,646080 5794900,646080 5797460))'
    polygon = ogr.CreateGeometryFromWkt(wkt)
    feature.SetGeometry(polygon)
    layer.CreateFeature(feature)
    feature = None
    data_source = None
           

4、shp拆分多个shp

import osgeo.ogr as ogr
import osgeo.osr as osr

def run():
    driver = ogr.GetDriverByName('ESRI Shapefile')
    fileName = "E:/cq/中小河流.shp"
    dataSource = driver.Open(fileName, 0)
    layer = dataSource.GetLayer(0)
    print("空间参考 :{0}".format(layer.GetSpatialRef()))
    for i in range(0, layer.GetFeatureCount()):
        feat = layer.GetFeature(i)
        wkt = feat.geometry()
        print(wkt)
        create_shp(i, str(wkt))


def create_shp(shp_name,wkt):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(f'E:/cq/tif/shp/{shp_name}.shp')
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    layer = data_source.CreateLayer(f'{shp_name}', srs, ogr.wkbMultiPolygon)
    feature = ogr.Feature(layer.GetLayerDefn())
    polygon = ogr.CreateGeometryFromWkt(wkt)
    feature.SetGeometry(polygon)
    layer.CreateFeature(feature)
    feature = None
    data_source = None
           

5、shp截取tif

from osgeo import gdal

    input_raster = r"E:/chqqh/tif/FA/o0001.tif"
    input_raster = gdal.Open(input_raster)
    input_shape = r"E:/chqqh/tif/shp/1.shp"  # or any other format
    output_raster = r'E:\chqqh\tif\test\o0001_1.tif'
    ds = gdal.Warp(output_raster,
                   input_raster,
                   format='GTiff',
                   cutlineDSName=input_shape,  # or any other file format
                   cutlineWhere="FIELD = 'whatever'",
                   # optionally you can filter your cutline (shapefile) based on attribute values
                   dstNodata=-9999)  # select the no data value you like
    ds = None