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 pdsf = 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 shapefilesf = 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 geometrylons, 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 polygonif 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 osrdef 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 = Nonedata_source = None

4、shp拆分多个shp

import osgeo.ogr as ogr
import osgeo.osr as osrdef 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 = Nonedata_source = None

5、shp截取tif

from osgeo import gdalinput_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 formatoutput_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 formatcutlineWhere="FIELD = 'whatever'",# optionally you can filter your cutline (shapefile) based on attribute valuesdstNodata=-9999)  # select the no data value you likeds = None

Python之shp文件相关推荐

C#调用python实现shp文件的CGCS2000坐标系向WGS84坐标系转换
前言:我们项目上最近有个需求,要把shp文件导入到我们的webgis系统中来,展示在前端的网页上.如果是任意的shp想展示在webgis中,我们就使用postgresql提供的ogr2ogr工具,直接 ...
Python打开shp文件报错：Set SHAPE_RESTORE_SHX config option to YES to restore or create it.
ERROR 4: Unable to open point.shx or point.SHX. Set SHAPE_RESTORE_SHX config option to YES to restor ...
python shp文件_对python 读取线的shp文件实例详解
如下所示: import shapefile sf = shapefile.reader("e:\\1.2\\cs\\dx_csl.shp") shapes = sf.shapes ...
shp文件批量定义投影坐标系python
这是官方示例,arcgis pro 对路径下的所有shp文件定义坐标系,转换坐标系类似(去上面的链接里获取). 代码注意事项: python的解释器要使用ArcGIS集成安装的,这样才能调用arcpy ...
Python|shp文件逐图斑自动生成地块界址点信息表——OGR库实现
工作中,有时需要从面shp文件中导出拐点坐标,并生成土地报批标准格式TXT文件或界址点信息表.ArcGIS中可以通过工具提取出拐点坐标,但是无法自动生成我们需要的成果格式,需要手动录入.今天分享基于P ...
python使用gdal将shp文件转为TIF
python使用gdal将shp文件转为TIF 方法一 # 缺少获取shp文件坐标系的步骤 def vector2raster(inputfilePath, outputfile, resp):sf ...
python处理shp和栅格文件的相关库shapefile、gdal等
读取含polygon的shp文件: def readshp(shp_path):sf = shapefile.Reader(shp_path)#创建reader类的对象进行shapefile文件的读取 ...
ncl 添加点shp文件_一：python读取shapefile文件
使用geopandas读取shapefile格式的矢量文件包括shapefile文件的读取.修改.保存 import geopandas as gpd # 导入包 segpath = r" ...
python使用gdal读取shp文件
from __future__ import division import pandas as pd from osgeo import ogr, osr, gdal import osdef re ...
Python 画降雨图，mash 数据掩膜裁剪出想要的区域，根据shp文件裁剪tif数据
1. 背景介绍 2. 导库 import numpy as np import xarray as xr import matplotlib.pyplot as plt from matplotlib ...