import rasterio
dataset = rasterio.open('forest_loss_porijogi_wgs84.tif')
print(dataset.name)
print(dataset.mode)
print(dataset.count)
print(dataset.width)
print(dataset.height)
print(dataset.crs)
print(dataset.bounds)
forest_loss_porijogi_wgs84.tif r 1 1326 687 EPSG:4326 BoundingBox(left=26.548502689, bottom=58.117685726, right=26.880002689, top=58.289435726)
print(dataset.profile)
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': 128.0, 'width': 1326, 'height': 687, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.0002500000000000014, 0.0, 26.548502689, 0.0, -0.0002500000000000043, 58.289435726), 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}
print(dataset.transform)
| 0.00, 0.00, 26.55| | 0.00,-0.00, 58.29| | 0.00, 0.00, 1.00|
for i in range(len(dataset.indexes) ):
print("{}: {}".format(i, dataset.dtypes[i]))
0: int16
# reading the first band (not from zero!)
band1 = dataset.read(1)
band1
array([[ 0, 0, 0, ..., 0, 0, 0], [ 0, 0, 0, ..., 0, 0, 0], [ 0, 0, 0, ..., 0, 0, 0], ..., [ 0, 17, 17, ..., 0, 0, 0], [ 0, 0, 0, ..., 0, 0, 0], [ 0, 0, 0, ..., 0, 0, 0]], dtype=int16)
import matplotlib.pyplot as plt
%matplotlib inline
plt.imshow(band1)
<matplotlib.image.AxesImage at 0x1e1992b02b0>
from rasterio.plot import show
show(dataset)
<matplotlib.axes._subplots.AxesSubplot at 0x1e1994ec630>
import numpy as np
# get classes
uniq_vals = np.unique(band1)
# display sorted order
print(sorted(uniq_vals))
# Patches = the matplotlib objects drawn
counts, bins = np.histogram(band1, bins=18)
# Print histogram outputs
for i in range(len(bins)-1):
print("bin lower bound:", bins[i])
print("counts:", counts[i])
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18] bin lower bound: 0.0 counts: 815440 bin lower bound: 1.0 counts: 4269 bin lower bound: 2.0 counts: 4179 bin lower bound: 3.0 counts: 1076 bin lower bound: 4.0 counts: 3273 bin lower bound: 5.0 counts: 4023 bin lower bound: 6.0 counts: 4418 bin lower bound: 7.0 counts: 3327 bin lower bound: 8.0 counts: 4611 bin lower bound: 9.0 counts: 3576 bin lower bound: 10.0 counts: 3983 bin lower bound: 11.0 counts: 5469 bin lower bound: 12.0 counts: 3706 bin lower bound: 13.0 counts: 2764 bin lower bound: 14.0 counts: 5237 bin lower bound: 15.0 counts: 6990 bin lower bound: 16.0 counts: 13465 bin lower bound: 17.0 counts: 21156
from matplotlib.patches import Patch
from matplotlib.colors import BoundaryNorm
from matplotlib import rcParams, cycler
fig, ax = plt.subplots()
cmap = plt.cm.viridis
lst = [int(x) for x in np.linspace(0,255,19)]
legend_patches = [Patch(color=icolor, label=label) for icolor, label in zip( cmap(lst), sorted(uniq_vals))]
ax.legend(handles=legend_patches,
facecolor="white",
edgecolor="white",
bbox_to_anchor=(1.35, 1))
plt.imshow(band1, cmap=cmap, interpolation='nearest')
<matplotlib.image.AxesImage at 0x1e19977e470>
from rasterio.plot import show_hist
show_hist(dataset, bins=19, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled', title="Histogram")
fig, (ax_dat, ax_hist) = plt.subplots(1, 2, figsize=(14,7))
ax_dat.legend(handles=legend_patches, facecolor="white", edgecolor="white")
show((dataset, 1), ax=ax_dat)
show_hist((dataset, 1), bins=19, ax=ax_hist)
plt.show()
dataset.close()
import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
dst_crs = 'EPSG:3301'
with rasterio.open('forest_loss_porijogi_wgs84.tif') as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open('forest_loss_porijogi_3301.tif', 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
with rasterio.open('forest_loss_porijogi_3301.tif', 'r') as data2:
show(data2, cmap=cmap)
display(data2.profile)
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': 128.0, 'width': 1107, 'height': 1088, 'count': 1, 'crs': CRS.from_epsg(3301), 'transform': Affine(18.286653237913317, 0.0, 649439.5163138957, 0.0, -18.286653237913317, 6464597.241364011), 'tiled': False, 'interleave': 'band'}
import geopandas as gpd
catchments = gpd.read_file('porijogi_sub_catchments.geojson')
display(catchments.crs)
display(catchments.head(5))
{'init': 'epsg:3301'}
OBJECTID | NAME_1 | AREA_1 | Shape_Leng | Shape_Area | ID | geometry | |
---|---|---|---|---|---|---|---|
0 | 8 | Idaoja | 3823.427995 | 35446.162219 | 3.823428e+07 | 1 | MULTIPOLYGON (((660834.858 6455555.914, 660851... |
1 | 9 | Keskjooks | 5087.809731 | 42814.174755 | 5.087810e+07 | 2 | MULTIPOLYGON (((666339.502 6455972.600, 666384... |
2 | 10 | Peeda | 5634.162684 | 47792.268153 | 5.634163e+07 | 3 | MULTIPOLYGON (((659914.002 6456514.131, 659817... |
3 | 11 | Sipe | 890.280919 | 16449.028656 | 8.902809e+06 | 4 | MULTIPOLYGON (((665928.914 6460634.243, 665985... |
4 | 12 | Tatra | 3306.643841 | 31108.960376 | 3.306644e+07 | 5 | MULTIPOLYGON (((658678.470 6457825.152, 658579... |
catchments.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x1e19beb0dd8>
import fiona
with fiona.open("porijogi_sub_catchments.geojson", "r") as vectorfile:
shapes = [feature["geometry"] for feature in vectorfile]
from rasterio.mask import mask
data2 = rasterio.open('forest_loss_porijogi_3301.tif', 'r')
# Clip the raster with Polygon
out_image, out_transform = mask(dataset=data2, shapes=shapes, crop=True)
out_meta = data2.meta.copy()
data2.close()
out_meta
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': 128.0, 'width': 1107, 'height': 1088, 'count': 1, 'crs': CRS.from_epsg(3301), 'transform': Affine(18.286653237913317, 0.0, 649439.5163138957, 0.0, -18.286653237913317, 6464597.241364011)}
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open("forest_loss_clipped.tif", "w", **out_meta) as dest:
dest.write(out_image)
with rasterio.open("forest_loss_clipped.tif", "r") as data3:
show(data3)
display(data3.profile)
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': 128.0, 'width': 1072, 'height': 1034, 'count': 1, 'crs': CRS.from_epsg(3301), 'transform': Affine(18.286653237913317, 0.0, 649896.6826448436, 0.0, -18.286653237913317, 6464158.3616863005), 'tiled': False, 'interleave': 'band'}
import geopandas as gpd
catchments = gpd.read_file('porijogi_sub_catchments.geojson')
demdata = rasterio.open('dem.tif')
print(demdata.name)
print(demdata.mode)
print(demdata.count)
print(demdata.width)
print(demdata.height)
print(demdata.crs)
print(demdata.bounds)
dem.tif r 1 1020 1065 EPSG:3301 BoundingBox(left=644405.5556, bottom=6443988.8889, right=675000.5556, top=6475948.8889)
fig, ax = plt.subplots(1, figsize=(9, 7))
show((demdata, 1), cmap='terrain', interpolation='none', ax=ax)
catchments.plot(ax=ax, facecolor="none", edgecolor='black', lw=0.7)
plt.title("Elevation in the Porijogi catchment")
plt.show()
from rasterstats import zonal_stats
zs = zonal_stats('porijogi_sub_catchments.geojson', 'dem.tif', stats=['mean','std'])
display(zs)
[{'mean': 108.81127854956439, 'std': 20.45245360418255}, {'mean': 86.88054631660887, 'std': 16.95209830197268}, {'mean': 122.27791004234241, 'std': 25.821034407200845}, {'mean': 76.41216968715197, 'std': 6.690340937594575}, {'mean': 82.89976042687574, 'std': 23.864544169477856}, {'mean': 66.31877151178183, 'std': 6.77909304322478}, {'mean': 112.26344569806881, 'std': 15.545179792343868}, {'mean': 59.404531339850486, 'std': 16.974837452873903}]
import pandas as pd
demstats_df = pd.DataFrame(zs)
demstats_df.rename(columns={'mean':'dem_mean','std':'dem_std'}, inplace=True)
catchments = pd.concat([catchments, demstats_df], axis=1)
fig, ax = plt.subplots(1, 1)
plt.title("Mean elevation per subcatchment")
catchments.plot(column='dem_mean', ax=ax, legend=True)
<matplotlib.axes._subplots.AxesSubplot at 0x1e19c1398d0>
import os
from glob import glob
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
# Stack the Landsat 8 bands
landsat_path = glob("Landsat8_pori__sr*.tif")
landsat_path.sort()
for idx, f in enumerate(landsat_path):
print(f"{idx}: {f}")
arr_st_l8, meta_l8 = es.stack(landsat_path, nodata=-9999)
0: Landsat8_pori__sr_band1.tif 1: Landsat8_pori__sr_band2.tif 2: Landsat8_pori__sr_band3.tif 3: Landsat8_pori__sr_band4.tif 4: Landsat8_pori__sr_band5.tif 5: Landsat8_pori__sr_band6.tif 6: Landsat8_pori__sr_band7.tif
b1 = arr_st_l8[4]
b2 = arr_st_l8[3]
ndvi_l8 = es.normalized_diff(b1, b2)
title = "Landsat 8 - Normalized Difference Vegetation Index (NDVI)"
# Turn off bytescale scaling due to float values for NDVI
ep.plot_bands(ndvi_l8, cmap="RdYlGn", cols=1, title=title, vmin=-1, vmax=1)
plt.tight_layout()
<Figure size 432x288 with 0 Axes>
# 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_l8, ndvi_class_bins)
# Apply the nodata mask to the newly classified NDVI data
ndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi_l8), ndvi_landsat_class)
# and check the number of unique values in our now classified dataset
np.unique(ndvi_landsat_class)
masked_array(data=[1, 2, 3, 4, 5], mask=False, fill_value=999999, dtype=int64)
# 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_l8 = np.unique(ndvi_landsat_class)
classes_l8 = classes_l8.tolist()
# The mask returns a value of none in the classes. remove that
classes_l8 = classes_l8[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_l8, 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()
# hillshading
import os
import numpy as np
import matplotlib.pyplot as plt
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
import rasterio as rio
dtm = 'dem.tif'
# Open the DEM with Rasterio
src = rio.open(dtm)
elevation = src.read(1)
# Plot the data
ep.plot_bands(elevation, cmap="gist_earth", title="DTM Without Hillshade", figsize=(10, 6) )
plt.tight_layout()
<Figure size 432x288 with 0 Axes>
hillshade = es.hillshade(elevation)
ep.plot_bands(hillshade, cbar=False, title="Hillshade made from DTM", figsize=(10, 6) )
plt.tight_layout()
<Figure size 432x288 with 0 Axes>
# Change the azimuth of the hillshade layer
hillshade_azimuth_210 = es.hillshade(elevation, azimuth=210)
# Plot the hillshade layer with the modified azimuth
ep.plot_bands(hillshade_azimuth_210, cbar=False, title="Hillshade with Azimuth set to 210 Degrees", figsize=(10, 6) )
plt.tight_layout()
<Figure size 432x288 with 0 Axes>
# Adjust the azimuth value
hillshade_angle_10 = es.hillshade(elevation, altitude=10)
# Plot the hillshade layer with the modified angle altitude
ep.plot_bands(hillshade_angle_10, cbar=False, title="Hillshade with Angle Altitude set to 10 Degrees", figsize=(10, 6) )
plt.tight_layout()
<Figure size 432x288 with 0 Axes>
# Plot the DEM and hillshade at the same time
fig, ax = plt.subplots(figsize=(10, 6))
ep.plot_bands(elevation, ax=ax, cmap="terrain", title="Lidar Digital Elevation Model (DEM)\n overlayed on top of a hillshade" )
ax.imshow(hillshade_angle_10, cmap="Greys", alpha=0.5)
plt.tight_layout()