Python矢量与栅格数据处理
原创 孟祥帅 祥帅的小屋 2022-10-21 18:55 发表于山东
持续更新ing，欢迎各位点赞关注哦！
1. Python处理栅格影像
1.1 读取tif
import rasterio
from rasterio.plot import show
from matplotlib import colors, cm
rs = 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 gpd
shp = 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 os
import geopandas
from shapely import geometry
import matplotlib.pyplot as plt
x1,y1=30,30
x2,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
Original: https://blog.csdn.net/qq_36980284/article/details/127452156
Author: 小孟的CDN
Title: Python与GIS
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/764380/
转载文章受原作者版权保护。转载请注明原作者出处！