-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathtask3.py
157 lines (130 loc) · 6.31 KB
/
task3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
from osgeo import gdal
import argparse
from task2 import handleTiff
import numpy as np
import ogr, os, osr
import matplotlib.pyplot as plt
import time
"""
Wider extent values:
minX = -2207050
maxX = -1002275
minY = -501171
maxY = -53369
"""
def readCommands():
"""
Read commandline arguments
"""
p = argparse.ArgumentParser(description=("Elevation change between two rasters"))
p.add_argument("--minX", dest ="minX", type=int, default=-2207050, help=("Minimum X bound"))
p.add_argument("--maxX", dest ="maxX", type=int, default=-1002275, help=("Maximum X bound"))
p.add_argument("--minY", dest ="minY", type=int, default=-501171, help=("Minimum Y Bound"))
p.add_argument("--maxY", dest ="maxY", type=int, default=-53369, help=("Maximum Y bound"))
cmdargs = p.parse_args()
return cmdargs
class clipTiff(object):
"""
Class to standardize rasters for further analysis
"""
def __init__(self,filename,minX,minY,maxX,maxY):
self.warpRaster(filename,minX,minY,maxX,maxY)
def warpRaster(self,filename,minX,minY,maxX,maxY,res=50):
"""
Function using gdal.Warp to reproject a raster
into a particular spatial extent for the EPSG:3031 projection
"""
ds = gdal.Open(str(filename))
self.out_filename = filename[:-4] + str("_clipped.tif") # creating output name
try:
print("--Standardizing--")
ds = gdal.Warp(str(self.out_filename), filename, xRes=res, yRes=-res, dstSRS='EPSG:3031', outputBounds=(minX,minY,maxX,maxY))
ds = None
except:
print("Sorry your bounds choices were invalid") # checking user chose valid bounds
return self.out_filename
class changeDetection(object):
"""
Class to detect change between two rasters
"""
def __init__(self,array1,array2):
self.arrayCalc(array1,array2) # load in array objects (from Tiffs) created in handleTiff class
self.volumnCalc(array1,array2)
def arrayCalc(self,array1,array2):
"""
Function to calculate the difference
between two raster arrays
"""
self.array_error = False # flag for calculus
if array1.data.shape != array2.data.shape:
print("You got an array problem boss") # arrays don't match each other
self.array_error = True # set flag
else:
self.array_difference = array2.data - array1.data # otherwise find their difference
self.reference = array2 = array1 # take each input array objects' attributes for this object
def volumnCalc(self,array1,array2):
"""
Function performing calculus on the arrays
to work out the total volumn of change
between them
"""
if self.array_error == False: # if the input arrays match in shape (can be used for calculus)
# get original resolution (pre-warping)
self.xRes = round(self.reference.pixelWidth)
self.yRes = -round(self.reference.pixelHeight)
# find the area of each array cell
self.cellsize = self.xRes*self.yRes # metres
# set up new array
self.vol_change_per_cell = np.copy(self.array_difference)
print("--Calculating volume change--")
# loop through extent
for i in np.arange(0,self.reference.nY-1):
for j in np.arange(0,self.reference.nX)-1:
if np.isnan(self.array_difference[i][j]) == False:
# calculate volumn change per cell
self.vol_change_per_cell[i][j] = self.array_difference[i][j] * self.cellsize
# total change (sum of all cells)
self.total_vol_change = np.nansum(self.vol_change_per_cell.flat)
# calculate total volumn per cell in each year
array2.vol_per_cell = array2.data * self.cellsize
array2.total_vol = np.nansum(array2.vol_per_cell.flat)
array1.vol_per_cell = array1.data * self.cellsize
array1.total_vol = np.nansum(array1.vol_per_cell.flat)
def writeTiff(self,array_to_write,filename):
"""
Take the filled elevation data and create an output raster
"""
# string addition to output filename
# set geolocation information (note geotiffs count down from top edge in Y)
geotransform = (self.reference.xOrigin, self.xRes, 0, self.reference.yOrigin, 0, -self.yRes)
# load data in to geotiff object
dst_ds = gdal.GetDriverByName('GTiff').Create(filename, self.reference.nX, self.reference.nY, 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference() # establish encoding
srs.ImportFromEPSG(3031) # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(array_to_write) # write image to the raster
dst_ds.GetRasterBand(1).SetNoDataValue(np.NaN) # set no data value
dst_ds.FlushCache() # write to disk
dst_ds = None
print("Image written to",filename)
if __name__=="__main__":
start_time = time.time()
cmd = readCommands()
# Two files being prepared for computation (standardisation prior to processing)
clip_file_1 = clipTiff(filename=r'./2009/2009_LVIS_dem_filled_200m.tif',minX=cmd.minX,minY=cmd.minY,maxX=cmd.maxX,maxY=cmd.maxY)
clip_file_2 = clipTiff(filename=r'./2015/2015_LVIS_dem_filled_200m.tif',minX=cmd.minX,minY=cmd.minY,maxX=cmd.maxX,maxY=cmd.maxY)
# Reading out the clipped raster arrays
array_1 = handleTiff(filename=str(clip_file_1.out_filename),readTiff=True)
array_2 = handleTiff(filename=str(clip_file_2.out_filename),readTiff=True)
# Reading these back in for calculations
output = changeDetection(array_1,array_2)
output.volumnCalc(array_1,array_2)
# Write out the raster of change detection
output.writeTiff(array_to_write=output.array_difference,filename=r'./results/elevation_change_output.tif')
print("--- %s seconds ---" % (time.time() - start_time))
plt.plot()
plt.ylabel("Change in Elevation (m)")
plt.xlabel("Relative Latitude (Top to Bottom)")
plt.title("Cross-sectional Change")
plt.savefig("xSection_Elev_Change.png")