تعریف شبکه جریان حوضه آبریز با پایتون و Pysheds + آموزش
اگر ما زمین پردازش با ماهیت GIS خود را به پایتون منتقل کنیم چه اتفاقی می افتد؟ چه اتفاقی می افتد اگر ما داده های فضایی رستری و برداری خود را به عنوان اشیاء و متغیرها در یک اسکریپت پایتون 3 پردازش کنیم؟ سپس می توانیم از خودمان بپرسیم که آیا لازم است که چرخه را مجددا بسازیم، لازم است یک جریان کاری را تغییر دهیم که قبلا در یک نرم افزار GIS کار می کند؟ دلیل و لزوم چیست؟
پاسخ ساده ای به این معضل وجود دارد: کنترل بیشتر
کار با Python به ما کنترل بیشتری بر روی geoprocessing خود می دهد، از آنجا که ما رابط کاربری گرافیکی (GUI) را با آیکون ها، دکمه ها و جعبه های محاوره ای آن را ترک می کنیم. با استفاده از پایتون در ژوپیتر نوت بوک، ما می توانیم با فایل های خاص پیوند داشته باشیم، geoprocess را تعریف کنیم و گزینه های آن، قطعه های داده های پیش نویس و نهایی، و نتایج خروجی به فرمت های SIG بردار / برش تولید کنیم. مزایای دیگری از تجزیه و تحلیل فضایی در پایتون وجود دارد؛ اینکه قابل بازیابی و سرعت پردازش هستند.
Pysheds پایتون 3 بسته طراحی شده برای تعریف حوضه و استخراج جریان شبکه است. در این کتابخانه نیاز به مجموعه ای از پردازش داده پیشرفته و کتابخانه های تجزیه و تحلیل فضایی به عنوان Numpy، Pandas، Scipy، Scikit-Image، Rasterio و... است. این آموزش روش کامل در یک نوت بوک Jupyter با پایتون و Pysheds را نشان می دهد:
- یک مدل ارتقاء دیجیتال (بدون غرقاب) وارد کنید
- تعیین شیب مسیر جریان
- تعریف حوضه آبریز
- تجزیه و تحلیل توده تجمعی جریان
- استخراج شبکه جریان
- تولید حوضه / بردار
از آنجا که اکثر کاربران در ویندوز کار می کنند، ما یک آموزش برای نصب Pyshed ها و کتابخانه های مورد نیاز در ویندوز با Anaconda انجام داده ایم.
لینک ها و دستورات مفید:
برای نصب Anaconda لطفا به اینجا مراجعه کنید.
صفحه اصلی Pysheds در اینجا در دسترس است.
برای نصب کتابخانه mplleaflet، لطفا در Anaconda Prompt تایپ کنید:
pip install mplleaflet
داده های ورودی:
شما می توانید داده های ورودی برای این آموزش را در این لینک دانلود کنید.
کد پایتون: (این کد را از اینجا دانلود کنید)
This is the python code used for the geoprocessing:
Import required libraries
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import geopandas as gpd
from pysheds.grid import Grid
import mplleaflet
%matplotlib inline
Open DEM file
grid = Grid.from_raster('../Rst/LocalDem.tif', data_name='dem')
def plotFigure(data, label, cmap='Blues'):
plt.figure(figsize=(12,10))
plt.imshow(data, extent=grid.extent, cmap=cmap)
plt.colorbar(label=label)
plt.grid()
plotFigure(grid.dem, 'Elevation (m)')
Define direction map
#N NE E SE S SW W NW
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
grid.flowdir(data='dem', out_name='dir', dirmap=dirmap)
plotFigure(grid.dir,'Flow Directiom','viridis')
Define catchment
# Specify pour point
x, y = 285612.017,2936416.682
# Delineate the catchment
grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
recursionlimit=15000, xytype='label', nodata_out=0)
# Clip the bounding box to the catchment
grid.clip_to('catch')
# Get a view of the catchment
demView = grid.view('dem', nodata=np.nan)
plotFigure(demView,'Elevation')
#export selected raster
grid.to_raster(demView, '../Output/clippedElevations.tif')
Define accumulation grid and stream network
grid.accumulation(data='catch', dirmap=dirmap, pad_inplace=False, out_name='acc')
accView = grid.view('acc', nodata=np.nan)
plotFigure(accView,"Cell Number",'PuRd')
streams = grid.extract_river_network('catch', 'acc', threshold=200, dirmap=dirmap)
streams["features"][:2]
def saveDict(dic,file):
f = open(file,'w')
f.write(str(dic))
f.close()
#save geojson as separate file
saveDict(streams,'../Output/streams.geojson')
Plot DEM and stream network
streamNet = gpd.read_file('../Output/streams.geojson')
streamNet.crs = {'init' :'epsg:32613'}
# The polygonize argument defaults to the grid mask when no arguments are supplied
shapes = grid.polygonize()
# Plot catchment boundaries
fig, ax = plt.subplots(figsize=(6.5, 6.5))
for shape in shapes:
coords = np.asarray(shape[0]['coordinates'][0])
ax.plot(coords[:,0], coords[:,1], color='cyan')
ax.set_xlim(grid.bbox[0], grid.bbox[2])
ax.set_ylim(grid.bbox[1], grid.bbox[3])
ax.set_title('Catchment boundary (vector)')
gpd.plotting.plot_dataframe(streamNet, None, cmap='Blues', ax=ax)
#ax = streamNet.plot()
mplleaflet.display(fig=ax.figure, crs=streamNet.crs, tiles='esri_aerial')
شناسه تلگرام مدیر سایت: SubBasin@
نشانی ایمیل: behzadsarhadi@gmail.com
(سوالات تخصصی را در گروه تلگرام ارسال کنید)
_______________________________________________________
نظرات (۰)