Python与GIS

Python矢量与栅格数据处理

原创 孟祥帅 祥帅的小屋 2022-10-21 18:55 发表于山东

持续更新ing，欢迎各位点赞关注哦！

python地理信息基础知识，基本图形的创建
https://github.com/CNFeffery/DataScienceStudyNotes

geopandas建立缓冲区与地图处理、合并、交叉、缓冲
https://blog.csdn.net/fengdu78/article/details/107054249

1. Python处理栅格影像

1.1 读取tif

import rasteriofrom rasterio.plot import showfrom matplotlib import colors, cmrs = rasterio.open(r'C:\Users\lenovo\Desktop\pzh_map_dispose\sf1.tif','r')result1=rs.read()

1.2 替换Nodata数据

#满足条件，则替换，否则保持原样result=np.where(result1==result1.min(),np.nan,result1)result# np.unique(rss[0])

out:

array([[[nan, nan, nan, ..., nan, nan, nan],

[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]]])

1.2 显示TIF

import geopandas as gpdshp = gpd.read_file(r"C:\Users\lenovo\Desktop\pzh_map_dispose\pzh_city.shp")fig, ax = plt.subplots(figsize=(5,9))shp.plot(ax=ax,color='none')show(result, transform=rs_mask.transform,ax=ax, cmap='gist_earth')fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.nanmin(result), vmax=np.nanmax(result)), cmap='gist_earth'), ax=ax,extend='both',fraction=0.05)

2. 从头开始的一个例子——完整版

2.1 创建shp面，并写入文件

import osimport geopandasfrom shapely import geometryimport matplotlib.pyplot as pltx1,y1=30,30x2,y2=50,50# 对应shapely.geometry中的Polygon，用于表示面，下面我们创建一个由若干Polygon对象组成cq = geopandas.GeoSeries([geometry.Polygon([(x1,y1), (x2,y1), (x2,y2), (x1,y2)]),geometry.Polygon([(x2,y1),(55,40), (x2,y2)])],index=['1', '2'], # 构建一个索引字段crs='EPSG:4326', # 坐标系是：WGS 1984)cq.to_file(r'simple_poly.shp',driver='ESRI Shapefile',encoding='utf-8')cq

2.2 使用geopandas读取

gdf=geopandas.read_file(r'simple_poly.shp')gdf

gdf.plot(column='index')

2.3 面转为栅格——根据字段转栅格

# 矢量数据转栅格:不需要模板def shpToRaster2(shp_path,raster_path,cellsize,field):'''矢量数据转栅格:不需要模板:param shp_path:需要转换的矢量数据:param raster_path:输出后栅格保存路径:param cellsize:栅格大小