Add remote sensing code

Finding-the-Nexus/Landsat8/EVI_calculation.py

@@ -0,0 +1,49 @@
+# -*- coding: utf-8 -*-
+"""
+Calculating Enhanced Vegetation Index (EVI)
+
+EVI is similar to Normalized Difference Vegetation Index (NDVI) and can be used to quantify vegetation greenness.
+However, EVI corrects for some atmospheric conditions and canopy background noise and is more sensitive in
+areas with dense vegetation. It incorporates an “L” value to adjust for canopy background,
+“C” values as coefficients for atmospheric resistance, and values from the blue band (B).
+These enhancements allow for index calculation as a ratio between the R and NIR values,
+while reducing the background noise, atmospheric noise, and saturation in most cases.
+
+EVI = G * ((NIR - R) / (NIR + C1 * R – C2 * B + L))
+
+- In Landsat 4-7, EVI = 2.5 * ((Band 4 – Band 3) / (Band 4 + 6 * Band 3 – 7.5 * Band 1 + 1)).
+
+- In Landsat 8, EVI = 2.5 * ((Band 5 – Band 4) / (Band 5 + 6 * Band 4 – 7.5 * Band 2 + 1)).
+
+Ref: https://www.usgs.gov/land-resources/nli/landsat/landsat-enhanced-vegetation-index?qt-science_support_page_related_con=0#qt-science_support_page_related_con
+
+"""
+
+from drought_prediction.utils.landsat8.landsat8_utils import read_band,image_show, image_histogram, image_adjust_brightness
+import numpy as np
+
+# """
+# Read bands
+# """
+b2 = read_band(2)
+b4 = read_band(4)
+b5 = read_band(5)
+b6 = read_band(6)
+
+
+# """
+# Calculate EVI
+# """
+
+# In Landsat 8, EVI = 2.5 * ((Band 5 – Band 4) / (Band 5 + 6 * Band 4 – 7.5 * Band 2 + 1)).
+
+# EVI = 2.5 * (NIR - RED ) / ( NIR + 6.0 * RED - 7.5 * BLUE+ 1.0 )
+
+EVI = 2.5 * (b5-b4)/(b5 + 6.0 * b4 - 7.5 * b2 + 1.0)
+
+img_ha = image_adjust_brightness(EVI, -0.7695, 1.459,
+                                 'RdYlGn',  #find other colourmaps @ https://matplotlib.org/tutorials/colors/colormaps.html
+                                 'Enhanced Vegetation Index (EVI)',
+                                 to_file='code-generated-EVI-rescaled.png')
+
+

Finding-the-Nexus/Landsat8/MNDWI_calculation.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+"""
+Calculating Modified Normalized Difference Water Index (MNDWI)
+
+MNDWI = (Green – SWIR) / (Green + SWIR)
+
+- For Landsat 7 data, NDWI = (Band 2 – Band 5) / (Band 2 + Band 5)
+
+- For Landsat 8 data, NDWI = (Band 3 – Band 6) / (Band 3 + Band 6)
+
+
+Ref: https://www.linkedin.com/pulse/ndvi-ndbi-ndwi-calculation-using-landsat-7-8-tek-bahadur-kshetri/
+
+"""
+
+from drought_prediction.utils.landsat8.landsat8_utils import read_band,image_show, image_histogram, image_adjust_brightness, get_gain_bias_angle
+
+# """
+# Read bands
+# """
+b3 = read_band(3)
+b4 = read_band(4)
+b5 = read_band(5)
+b6 = read_band(6)
+
+
+# """
+# Calculate MNDWI
+# """
+
+
+# For Landsat 8 data, MNDWI = (Band 3 – Band 6) / (Band 3 + Band 6)
+
+MNDWI = (b3 - b6)/(b3 + b6)
+
+img_ha = image_adjust_brightness(MNDWI,
+                                 -0.7695, 1.459,
+                                 'Blues',  #find other colourmaps https://matplotlib.org/tutorials/colors/colormaps.html
+                                 'Modified Normalized Difference Water Index (MNDWI)',
+                                 to_file='code-generated-MNDWI.png')
+
+

Finding-the-Nexus/Landsat8/NDVI_calculation.py

@@ -0,0 +1,123 @@
+# -*- coding: utf-8 -*-
+"""
+Calculating Normalized Difference Vegetation Index (NDVI)
+
+The normalized difference vegetation index (NDVI) is used to assess the state of vegetation.
+In living plants chlorophyll-A, from the photosynthetic machinery, strongly absorbs red color;
+on the other hand, near-infrared light is strongly reflected.
+Live, healthy vegetation reflects around 8% of red light and 50% of near-infrared light.
+Dead, unhealthy, or sparse vegetation reflects approximately 30% of red light and 40% of near-infrared light.
+
+By its formulation, NDVI ranges from -1 to +1. In practice, an area of an image containing living vegetation will have NDVI in the range 0.3 to 0.8.
+High water content clouds and snow will have negative values of the index. Bodies of water, having low reflectance in both Band 4 and 5,
+exhibit very low positive or negative index. Soil, having slightly higher reflectance in near-infrared than in red, will produce low positive values of the index.
+
+- In Landsat 8, NDVI = (Band 5 – Band 4) / (Band 5 + Band 4).
+
+Ref: https://nbviewer.jupyter.org/github/HyperionAnalytics/PyDataNYC2014/blob/master/ndvi_calculation.ipynb
+
+"""
+from drought_prediction.utils.landsat8.landsat8_utils import read_band,image_show, image_histogram, image_adjust_brightness, get_gain_bias_angle
+import numpy as np
+from matplotlib.colors import ListedColormap
+import matplotlib.pyplot as plt
+import earthpy.plot as ep
+
+
+
+# """
+# Read bands
+# """
+b4 = read_band(4)
+b5 = read_band(5)
+
+# """
+# Get reflectance gain and bias, and Sun elevation angle; calculate top-of-atmosphere reflectance
+# """
+#
+#
+# b4_gain, b4_bias, angle = get_gain_bias_angle(4)
+# b5_gain, b5_bias, angle = get_gain_bias_angle(5)
+#
+# b4_lambda_refl  = (b4_gain * b4 + b4_bias) / np.sin(angle)
+# b5_lambda_refl  = (b5_gain * b5 + b5_bias) / np.sin(angle)
+#
+#
+# """
+# Calculate NDVI
+# """
+#
+# ndvi = (b5_lambda_refl - b4_lambda_refl) / (b5_lambda_refl + b4_lambda_refl)
+# img_ha = image_adjust_brightness(ndvi, -0.7695, 1.459, 'OrRd', 'NDVI')
+
+# In Landsat 8, NDVI = (Band 5 – Band 4) / (Band 5 + Band 4).
+
+ndvi = (b5-b4)/(b5+b4)
+
+print('NDVI min value: ',ndvi.min())
+print('NDVI max value: ',ndvi.max())
+
+
+img_ha = image_adjust_brightness(ndvi,
+                                 -0.7695, 1.459,
+                                 'RdYlGn',  #find other colourmaps https://matplotlib.org/tutorials/colors/colormaps.html
+                                 'Normalized Difference Vegetation Index (NDVI)',
+                                 to_file='code-generated-NDVI2.png')
+
+
+
+# Reference for the following: https://earthpy.readthedocs.io/en/latest/gallery_vignettes/plot_calculate_classify_ndvi.html
+
+"""
+Classify NDVI - Categorise (or classify) the NDVI results into useful classes. 
+Values under 0 will be classified together as no vegetation. 
+Additional classes will be created for bare area and low, moderate, and high vegetation areas.
+"""
+
+# Create classes and apply to NDVI results
+ndvi_class_bins = [-np.inf, 0, 0.1, 0.25, 0.4, np.inf]
+ndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)
+
+# Apply the nodata mask to the newly classified NDVI data
+ndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi), ndvi_landsat_class)
+np.unique(ndvi_landsat_class)
+
+
+"""
+Plot Classified NDVI With Categorical Legend - EarthPy Draw_Legend()
+"""
+
+# Define color map
+nbr_colors = ["gray", "y", "yellowgreen", "g", "darkgreen"]
+nbr_cmap = ListedColormap(nbr_colors)
+
+# Define class names
+ndvi_cat_names = [
+    "No Vegetation",
+    "Bare Area",
+    "Low Vegetation",
+    "Moderate Vegetation",
+    "High Vegetation",
+]
+
+# Get list of classes
+classes = np.unique(ndvi_landsat_class)
+classes = classes.tolist()
+# The mask returns a value of none in the classes. remove that
+classes = classes[0:5]
+
+# Plot your data
+fig, ax = plt.subplots(figsize=(12, 12))
+im = ax.imshow(ndvi_landsat_class, cmap=nbr_cmap)
+
+ep.draw_legend(im_ax=im, classes=classes, titles=ndvi_cat_names)
+ax.set_title(
+    "Landsat 8 - Normalized Difference Vegetation Index (NDVI) Classes",
+    fontsize=14,
+)
+ax.set_axis_off()
+
+# Auto adjust subplot to fit figure size
+plt.tight_layout()
+
+plt.show()

Finding-the-Nexus/Landsat8/NDVI_change_over_time.py

@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+"""
+Calculate change in NDVI over time
+
+- In Landsat 8, NDVI = (Band 5 – Band 4) / (Band 5 + Band 4).
+
+Ref: https://geohackweek.github.io/raster/04-workingwithrasters/
+
+"""
+from drought_prediction.utils.landsat8.landsat8_utils import read_band,image_show, image_histogram, image_adjust_brightness, get_gain_bias_angle
+
+
+
+# """
+# Read bands
+# """
+b4 = read_band(4)
+b5 = read_band(5)
+
+# """
+# Get reflectance gain and bias, and Sun elevation angle; calculate top-of-atmosphere reflectance
+# """
+#
+#
+# b4_gain, b4_bias, angle = get_gain_bias_angle(4)
+# b5_gain, b5_bias, angle = get_gain_bias_angle(5)
+#
+# b4_lambda_refl  = (b4_gain * b4 + b4_bias) / np.sin(angle)
+# b5_lambda_refl  = (b5_gain * b5 + b5_bias) / np.sin(angle)
+#
+#
+# """
+# Calculate NDVI
+# """
+#
+# ndvi = (b5_lambda_refl - b4_lambda_refl) / (b5_lambda_refl + b4_lambda_refl)
+# img_ha = image_adjust_brightness(ndvi, -0.7695, 1.459, 'OrRd', 'NDVI')
+
+# In Landsat 8, NDVI = (Band 5 – Band 4) / (Band 5 + Band 4).
+
+ndvi = (b5-b4)/(b5+b4)
+img_ha = image_adjust_brightness(ndvi,
+                                 -0.7695, 1.459,
+                                 'RdYlGn',  #find other colourmaps https://matplotlib.org/tutorials/colors/colormaps.html
+                                 'Normalized Difference Vegetation Index (NDVI)',
+                                 to_file='code-generated-NDVI2.png')
+
+

Finding-the-Nexus/Landsat8/NDVI_stacked_version.py

@@ -0,0 +1,96 @@
+from matplotlib import pyplot as plt
+import earthpy.spatial as es
+import rasterio as rio
+import numpy as np
+from matplotlib.colors import ListedColormap
+import matplotlib.pyplot as plt
+import earthpy.plot as ep
+
+from Landsat8.landsat8_utils import get_list_of_data
+
+land_stack, \
+land_meta, \
+landsat_post_process_path = get_list_of_data('/Users/gkalliatakis/Desktop/Gedo/LC081650582018022501T1-SC20200416194144')
+
+
+# Once you have stacked your data, you can import it and work with it as you need to!
+with rio.open(landsat_post_process_path) as src:
+    landsat_stacked_data = src.read(masked=True)
+    csf_meta = src.meta
+
+
+nir_band = landsat_stacked_data[5]
+red_band = landsat_stacked_data[4]
+
+# Calculate normalized difference vegetation index
+ndvi = es.normalized_diff(b1=nir_band, b2=red_band)
+
+
+titles = ["Landsat 8 - Normalized Difference Vegetation Index (NDVI)"]
+
+# Turn off bytescale scaling due to float values for NDVI
+ep.plot_bands(ndvi,
+              cmap="RdYlGn",
+              cols=1,
+              title=titles,
+              scale=False,
+              vmin=-1, vmax=1
+)
+
+
+
+
+# Reference for the following: https://earthpy.readthedocs.io/en/latest/gallery_vignettes/plot_calculate_classify_ndvi.html
+
+"""
+Classify NDVI - Categorise (or classify) the NDVI results into useful classes. 
+Values under 0 will be classified together as no vegetation. 
+Additional classes will be created for bare area and low, moderate, and high vegetation areas.
+"""
+
+# Create classes and apply to NDVI results
+ndvi_class_bins = [-np.inf, 0, 0.1, 0.25, 0.4, np.inf]
+ndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)
+
+# Apply the nodata mask to the newly classified NDVI data
+ndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi), ndvi_landsat_class)
+np.unique(ndvi_landsat_class)
+
+
+"""
+Plot Classified NDVI With Categorical Legend - EarthPy Draw_Legend()
+"""
+
+# Define color map
+nbr_colors = ["gray", "y", "yellowgreen", "g", "darkgreen"]
+nbr_cmap = ListedColormap(nbr_colors)
+
+# Define class names
+ndvi_cat_names = [
+    "No Vegetation",
+    "Bare Area",
+    "Low Vegetation",
+    "Moderate Vegetation",
+    "High Vegetation",
+]
+
+# Get list of classes
+classes = np.unique(ndvi_landsat_class)
+classes = classes.tolist()
+# The mask returns a value of none in the classes. remove that
+classes = classes[0:5]
+
+# Plot your data
+fig, ax = plt.subplots(figsize=(12, 12))
+im = ax.imshow(ndvi_landsat_class, cmap=nbr_cmap)
+
+ep.draw_legend(im_ax=im, classes=classes, titles=ndvi_cat_names)
+ax.set_title("Landsat 8 - Normalized Difference Vegetation Index (NDVI) Classes",
+             fontsize=14,
+)
+ax.set_axis_off()
+
+# Auto adjust subplot to fit figure size
+plt.tight_layout()
+
+plt.show()