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