تجزیه و تحلیل تغییر پوشش زمین با پایتون و 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
داده ها ورودی بالا را از اینجا دانلود کنید.
تماس واتساپ: 09190622992 (98+)
شناسه تلگرام مدیر سایت: Basin_Ir_bot@
نشانی ایمیل: www.Basin.ir@Gmail.com
(سوالات تخصصی را در گروه تلگرام ارسال کنید)
_______________________________________________________
تماس با ما
تماس واتساپ
گروه تلگرام ما
هشدارهای ما
آمار و اطلاعات
سفارش پروژه
در منابع آب
نظرات (۰)