-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathLab05_2.4.3.py
151 lines (106 loc) · 5.35 KB
/
Lab05_2.4.3.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
############################################################
# 2.4.3 Raster Tutorials / Fire Analysis
#
#Provides functions to:
# 2.4.3 Identify High/Med/Low Heat Thresholds of Mills Fire
# > Raster Information
# > Clip To Specified Bounds
# > Downsample
# > Descriminate
# > Combine
#
#Information:
# Fire intensity models have utility for various industries.
# Cal Fire Intel uses these inputs to generate predictive
# models that assist in the determination of where the fire
# will grow, and how fast. Post fire damage assesment teams
# may use this data to plan rehabilitation of devistated
# areas based on burn intensity upon exposed flora.
#
# By: Tony Ramos
# Date: 02/19/2023
#
############################################################
import sys
# Open source spatial libraries
import shapely
import numpy as np
from osgeo import gdal
import math
import random
import os
# SpaPy libraries
from SpaPy import SpaBase
from SpaPy import SpaPlot
from SpaPy import SpaVectors
from SpaPy import SpaView
from SpaPy import SpaReferencing
from SpaPy import SpaDensify
from SpaPy import SpaView
from SpaPy import SpaRasters
from SpaPy import SpaTopo
from SpaPy import SpaRasterVectors
######################################################################
# File Paths
# A DEM of Mt St Helens before the eruption
HotspotHighlight="SpaPyTests/Data/MillsFire/LWIR_QuickMosaic_Hotspot_Highlight_9327_4x.tiff"
# A DEM of Mt St Helens after the eruption
LWIRFullRez="SpaPyTests/Data/MillsFire/LWIR_QuickMosaic_16-bit_9327.tiff"
# A temporary folder for outputs
TempFolderPath3="SpaPyTests/Temp3/"
os.makedirs(TempFolderPath3)
######################################################################
# RASTER INFORMATION
# Define TheDataset with our SpaDatasetRaster function
LWIR_Raster=SpaRasters.SpaDatasetRaster()
# Load our raster file in as TheDataset
LWIR_Raster.Load(LWIRFullRez)
# Get the dimensions of the raster in pixels
print("Width in pixels: "+format(LWIR_Raster.GetWidthInPixels()))
print("Height in pixels: "+format(LWIR_Raster.GetHeightInPixels()))
# Get the number of bands of data in the raster (e.g. DEMs have 1, RGB data has 3)
print("Number of Bands: "+format(LWIR_Raster.GetNumBands()))
# Get the GDAL type of pixels which includes GDT_Int16, GDT_Int32, GDT_Float32 and GDT_Float64
print("Pixel Type: "+format(LWIR_Raster.GetType()))
# Get the spatial reference
print("Coordinate Reference System/Spatial Reference: "+format(LWIR_Raster.GetCRS()))
# Get the resolution or dimensions of each pixel in the raster
print("Resolution (x,y): "+format(LWIR_Raster.GetResolution()))
# Get the spatial bounds of the raster
TheBounds=LWIR_Raster.GetBounds()
print("TheBounds (XMin,YMin,XMax,YMax): "+format(TheBounds))
# Get additional information on the raster
TheBandStats=LWIR_Raster.GetBandInfo(1)
print("TheBandStats="+format(TheBandStats))
# Get one band of data from the raster
TheBand=LWIR_Raster.GetBand(0)
print("TheBands: "+format(TheBand))
SpaView.Show(LWIR_Raster) # View original raster
#####################################################################
# CLIP TO SPECIFIED BOUNDS
# Use SpaRaster Crop tool and input raster for clip and desired bounds
ClippedRaster=SpaRasters.Crop(LWIR_Raster,[-122.415767,41.40831,-122.37479,41.5014]) # Bounds set for SpaView. Actual bounds dont view properly and are set later.
ClippedRaster.Save(TempFolderPath3+"Cropped.tif") # Save Result
SpaView.Show(ClippedRaster) # View Fire Area
#####################################################################
# DOWNSAMPLE
# Resample (Resample variable is the denominator, Original resolution is the numerator. Resolution is devided by the input value)
DownSample=SpaRasters.Resample(ClippedRaster,0.28) # Downsample raster to 5 meter spatial resolution.
DownSample.Save(TempFolderPath3 + "DownSampledRaster.tif") # Save Result
#####################################################################
# DESCRIMINATE
# select for high, medium, low heat thresholds.
LowHeat=SpaRasters.GreaterThanOrEqual(TempFolderPath3 + "DownSampledRaster.tif",33000) # Values below 33000 16-Bit Radiometric Resolution are no fire (0).
MediumHeat=SpaRasters.GreaterThanOrEqual(TempFolderPath3 + "DownSampledRaster.tif",39000) # Values between 39000-53000 are medium heat.
HighHeat=SpaRasters.GreaterThanOrEqual(TempFolderPath3 + "DownSampledRaster.tif",53000) # Values between 53000-65536(Max) are high heat
#####################################################################
# COMBINE
# Raster math to combine descriminated values.
FireClass = HighHeat + MediumHeat + LowHeat # Sum of all heats where 0 = NoHeat; 0.001-1 = LowHeat; 1.001-2 = MediumHeat; 2.001-3 = HighHeat
FireClass.Save(TempFolderPath3 + "Fire_Class.tif") # Save Result
FireClassCrop = SpaRasters.Crop(TempFolderPath3 + "Fire_Class.tif",[-122.415767,41.43,-122.37479,41.5014]) # Re-Cropped for actual bound
FireClassCrop.Save(TempFolderPath3 + "Fire_Class_Final.tif") # Saved Final Output
SpaView.Show(FireClass) # View heat classes.
SpaView.Show(FireClassCrop) # View to illustrate that corrected bounds doesn't render right in SpaView. Output raster bounds are correct.
######################################################################
# END