تجزیه و تحلیل تغییر پوشش زمین با پایتون و GDAL :: بیسین - سایت تخصصی مهندسی آب

ابزار وبمستر

Bootstrap Example

عضويت در خبرنامه ايـميـل پايگاه بيسيــن - عضويت پس از کليک بر روي لينک فعال سازي که براي شما ارسال خواهد شد تکميل مي شود

پشتيباني شده با بيسين

تجزیه و تحلیل تغییر پوشش زمین با پایتون و GDAL


تصاویر ماهواره ای توانایی دیدن سطح زمین در سال های اخیر را برای ما فراهم کرده است، اما ما در درک پویایی پوشش زمین و تعامل با عوامل اقتصادی، جامعه شناختی و سیاسی چندان موفق نبوده ایم. برخی از نقص های این تجزیه و تحلیل در استفاده از نرم افزار تجاری GIS مشاهده شد، اما محدودیت های دیگری نیز در نحوه اعمال فرآیندهای منطقی و ریاضی بر روی مجموعه ای از تصاویر ماهواره ای وجود دارد. کار با داده های جغرافیایی روی پایتون امکان فیلتر کردن، محاسبه، برش، تلفیق و صادرات داده های رستری و برداری را با استفاده شرایط کارآمدی از توان محاسباتی فراهم می کند و دامنه بیشتری در تجزیه و تحلیل داده ها دارد.


این آموزش روش کامل ایجاد یک لایه تغییر سطح پوشش زمین از مقایسه شاخص های پوشش گیاهی تولید شده (NDVI) با استفاده از کتابخانه های پایتون و Numpy و GDAL را نشان می دهد. محوطه تغییر پوشش زمین که در آنجا با برخی از ابزارهای GDAL و Osgeo تولید شده و تجزیه و تحلیل جنگل زدایی بر اساس داده های خروجی و تصاویر تاریخی از Google Earth انجام شده است.


کد پایتون

این کد پایتون برای آموزش است:


واردات به کتابخانه نیاز دارد

در این آموزش فقط از کتابخانه های اصلی Python ،Scipy ،Numpy ،Matplotlib و غیره و همچنین GDAL استفاده شده است. برای کاربران ویندوز، موثرترین روش بارگیری چرخ GDAL از اینجا و نصب از طریق روش pip است.


from osgeo import ogr, gdal, osr

import numpy as np

import os

import matplotlib.pyplot as plt


تنظیم پرونده های ورودی و خروجی

ما مسیر باند های رستری و NDVI خروجی و پوشش های تغییر سطح زمین را اعلام می کنیم. مسیر شکل تغییر پوشش زمین نیز مشخص شده است.


#Input Raster and Vector Paths


#Image-2019

path_B5_2019="../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B5_clip.TIF"

path_B4_2019="../Image20190203clip/LC08_L1TP_029047_20190203_20190206_01_T1_B4_clip.TIF"


#Image-2014

path_B5_2014="../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B5_clip.TIF"

path_B4_2014="../Image20140205clip/LC08_L1TP_029047_20140205_20170307_01_T1_B4_clip.TIF"


#Output Files


#Output NDVI Rasters 

path_NDVI_2019 = '../Output/NDVI2019.tif'

path_NDVI_2014 = '../Output/NDVI2014.tif'

path_NDVIChange_19_14 = '../Output/NDVIChange_19_14.tif'


#NDVI Contours

contours_NDVIChange_19_14 = '../Output/NDVIChange_19_14.shp'


باندهای تصویر Landsat را با GDAL باز کنید

در این قسمت قرمز (باند 4) و نزدیک NIR مادون قرمز (باند 5) را با دستورات کتابخانه GDAL باز می کنیم و سپس تصاویر را به صورت آرایه های ماتریس با نوع Float 32 بیت می خوانیم.



#Open raster bands

B5_2019 = gdal.Open(path_B5_2019)

B4_2019 = gdal.Open(path_B4_2019)

B5_2014 = gdal.Open(path_B5_2014)

B4_2014 = gdal.Open(path_B4_2014)


#Read bands as matrix arrays

B52019_Data = B5_2019.GetRasterBand(1).ReadAsArray().astype(np.float32)

B42019_Data = B4_2019.GetRasterBand(1).ReadAsArray().astype(np.float32)

B52014_Data = B5_2014.GetRasterBand(1).ReadAsArray().astype(np.float32)

B42014_Data = B4_2014.GetRasterBand(1).ReadAsArray().astype(np.float32)


اندازه های ماتریس و پارامترهای دگرسانس را با یکدیگر مقایسه کنید

عملکرد بین باند Landsat 8 بستگی به این دارد که پارامترهای وضوح، اندازه، src و geotransformation برای هر دو تصویر یکسان باشد. در صورتی که این ویژگی ها با پیچ و خمی یکسان نباشند، نیاز به دفع مجدد، مقیاس یا هر یک از فرآیند های جغرافیایی دیگر ضروری است.


print(B5_2014.GetProjection()[:80])

print(B5_2019.GetProjection()[:80])

if B5_2014.GetProjection()[:80]==B5_2019.GetProjection()[:80]: print('SRC OK')


PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84

PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84

SRC OK


print(B52014_Data.shape)

print(B52019_Data.shape)

if B52014_Data.shape==B52019_Data.shape: print('Array Size OK')


(610, 597)

(610, 597)

Array Size OK


print(B5_2014.GetGeoTransform())

print(B5_2019.GetGeoTransform())

if B5_2014.GetGeoTransform()==B5_2019.GetGeoTransform(): print('Geotransformation OK')


(652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0)

(652500.0, 29.98324958123953, 0.0, 2166000.0, 0.0, -30.0)

Geotransformation OK


پارامترهای جغرافیایی را دریافت کنید

از آنجا که ما مقایسه کرده ایم که گروه های رستر اندازه آرایه یکسانی دارند و در یک موقعیت قرار دارند می توانیم برخی پارامترهای مکانی را بدست آوریم.


geotransform = B5_2014.GetGeoTransform()


originX,pixelWidth,empty,finalY,empty2,pixelHeight=geotransform

cols =  B5_2014.RasterXSize

rows =  B5_2014.RasterYSize


projection = B5_2014.GetProjection()


finalX = originX + pixelWidth * cols

originY = finalY + pixelHeight * rows


NDVI را محاسبه کرده و آنها را به عنوان پرونده ای دیگر ذخیره کنید

ما می توانیم NDVI را به عنوان یک جبر ماتریس با برخی از توابع Numpy محاسبه کنیم. با استفاده از GDAL و پارامترهای دگرگونی می توان نتیجه را به عنوان یک خروجی صادر کرد. برای داشتن کد مؤثرتر و واضح تر، ما یک تابع برای صادرات شناسه ها ایجاد می کنیم.


ndvi2014 = np.divide(B52014_Data - B42014_Data, B52014_Data+ B42014_Data,where=(B52014_Data - B42014_Data)!=0)

ndvi2014[ndvi2014 == 0] = -999


ndvi2019 = np.divide(B52019_Data - B42019_Data, B52019_Data+ B42019_Data,where=(B52019_Data - B42019_Data)!=0)

ndvi2019[ndvi2019 == 0] = -999


def saveRaster(dataset,datasetPath,cols,rows,projection):

    rasterSet = gdal.GetDriverByName('GTiff').Create(datasetPath, cols, rows,1,gdal.GDT_Float32)

    rasterSet.SetProjection(projection)

    rasterSet.SetGeoTransform(geotransform)

    rasterSet.GetRasterBand(1).WriteArray(dataset)

    rasterSet.GetRasterBand(1).SetNoDataValue(-999)

    rasterSet = None


saveRaster(ndvi2014,path_NDVI_2014,cols,rows,projection)


saveRaster(ndvi2019,path_NDVI_2019,cols,rows,projection)


تصاویر طرح NDVI

ما یک تابع برای ترسیم تصاویر NDVI حاصل با یک colorbar ایجاد می کنیم


extentArray = [originX,finalX,originY,finalY]

def plotNDVI(ndviImage,extentArray,vmin,cmap):

    ndvi = gdal.Open(ndviImage)

    ds2019 = ndvi.ReadAsArray()

    plt.figure(figsize=(20,15))

    im = plt.imshow(ds2019, vmin=vmin, cmap=cmap, extent=extentArray)#

    plt.colorbar(im, fraction=0.015)

    plt.xlabel('Este')

    plt.ylabel('Norte')

    plt.show()


plotNDVI(path_NDVI_2014,extentArray,0,'YlGn')



plotNDVI(path_NDVI_2019,extentArray,0,'YlGn')



یک تصویر تغییر پوشش زمین ایجاد کنید

ما براساس تفاوت های NDVI از تصاویر سال 2019 و 2014 یک تصویر تغییر پوشش زمین ایجاد می کنیم. تصویر به صورت یک پرونده جداگانه ذخیره می شود.


ndviChange = ndvi2019-ndvi2014

ndviChange = np.where((ndvi2014>-999) & (ndvi2019>-999),ndviChange,-999)

ndviChange


array([[-9.9900000e+02,  3.5470784e-02,  3.7951261e-02, ...,

         3.1590372e-02,  3.8002759e-02,  2.6692629e-02],

       [-9.9900000e+02,  3.7654012e-02,  5.8923483e-02, ...,

        -5.8691800e-03,  1.8922418e-02,  1.9506305e-02],

       [-9.9900000e+02,  3.2184020e-02,  3.7395060e-02, ...,

        -6.3773081e-02, -3.0675709e-02,  3.8942695e-02],

       ...,

       [-9.9900000e+02, -1.6597062e-02,  1.1402190e-02, ...,

         2.7218461e-03,  2.5526762e-02,  6.6639006e-02],

       [-9.9900000e+02,  5.9944689e-03, -4.0770471e-03, ...,

         5.7501763e-02,  4.6757758e-03,  8.4484324e-02],

       [-9.9900000e+02, -9.9900000e+02, -9.9900000e+02, ...,

        -9.9900000e+02, -9.9900000e+02, -9.9900000e+02]], dtype=float32)


saveRaster(ndviChange,path_NDVIChange_19_14,cols,rows,projection)


plotNDVI(path_NDVIChange_19_14,extentArray,-0.5,'Spectral')



خطوط کانتور را ایجاد کنید

سرانجام، خطوط کانتور از مقادیر NDVI تولید می شوند. این کار با ابزار کتابخانه GDAL انجام می شود.


Dataset_ndvi = gdal.Open(path_NDVIChange_19_14)#path_NDVI_2014

ndvi_raster = Dataset_ndvi.GetRasterBand(1)


ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(contours_NDVIChange_19_14)


prj=Dataset_ndvi.GetProjectionRef()#GetProjection()


srs = osr.SpatialReference(wkt=prj)#

#srs.ImportFromProj4(prj)


contour_shp = ogr_ds.CreateLayer('contour', srs)

field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)

contour_shp.CreateField(field_defn)

field_defn = ogr.FieldDefn("ndviChange", ogr.OFTReal)

contour_shp.CreateField(field_defn)

#Generate Contourlines

gdal.ContourGenerate(ndvi_raster, 0.1, 0, [], 1, -999, contour_shp, 0, 1)

ogr_ds = None


داده ها ورودی بالا را از اینجا دانلود کنید.



×

راهنماي حل مشکل دانلود: با توجه به مسدود شدن درايو گوگل در ايران از آي پي ديگر کشورها براي دانلود فايل ها استفاده کنيد.




سفارش پروژه داريد؟ يا قصد همکاري در انجام پروژه؟ و يا قصد فروش فايل خود؟

با فشردن دکمه زير يکي از بخش هاي "سفارش انجام پروژه" يا "همکاري با بيسين" و يا "فروش فايل" را انتخاب فرماييد







آمار آنلاين-مقايسه اي بارش در حوضه هاي اصلي کشور

منبع: وزارت نيرو - اين نمودار ممکن است براي دقايقي به دليل بروزرساني غير فعال شود


آخرين تصوير ماهواره هواشناسي - موقعيت ايران

W3Schools


نظرات  (۰)

فرم ارسال نظر

ارسال نظر آزاد است، اما اگر قبلا در بیان ثبت نام کرده اید می توانید ابتدا وارد شوید.
شما میتوانید از این تگهای html استفاده کنید:
<b> یا <strong>، <em> یا <i>، <u>، <strike> یا <s>، <sup>، <sub>، <blockquote>، <code>، <pre>، <hr>، <br>، <p>، <a href="" title="">، <span style="">، <div align="">

درباره بهترين هاي بيسيـــن بدانيد...

Bird

يکي از مهمترين اهداف اين سايت تهيه آموزش هاي روان از ابزارهاي کاربردي علوم آب است.

اهميت مطالعات محيطي با ابزارهاي نوين در چيست؟

امروز با فارغ التحصيلي جمع کثير دانشجويان سالهاي گذشته و حال، با گذر از کمي گرايي ديگر صرف وجود مدارک دانشگاهي حرف اول را در بازار کار نمي زند؛ بلکه سنجش ديگري ملاک؛ و شايسته سالاري به ناچار! باب خواهد شد. يکي از مهم ترين لوازم توسعه علمي در هر کشور و ارائه موضوعات ابتکاري، بهره گيري از ابزار نوين است، بيسين با همکاري مخاطبان مي تواند در حيطه علوم آب به معرفي اين مهم بپردازد.

جستجو در بيسين


ابزارهاي نوين

بیسین - سایت تخصصی مهندسی آب

بیسین جهت ارائه مطالب و خدمات تخصصی در حیطه نرم افزارها و مدل های شبیه سازی مهندسی آب با رویکرد پژوهشی-آموزشی ایجاد شده است که توسعه خود را در گرو همکاری مخاطبان می بیند.

اطلاعات سايت

  • www.Basin.ir@gmail.com
  • بهزاد سرهادي
  • تاريخ امروز:
  • شناسه تلگرام: SubBasin
  • شماره واتساپ: 09190622992-098
  • شماره تماس: 09190622992-098

W3Schools