-
Notifications
You must be signed in to change notification settings - Fork 3
/
show_raster.py
82 lines (59 loc) · 2.58 KB
/
show_raster.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
#!/usr/bin/env python
""" show_plot.py
Plot raster data of the CSA MET for demonstration/visualization purposes
Author: Olivier Lamarre <[email protected]>
Affl.: Space and Terrestrial Autonomous Robotic Systems Laboratory
University of Toronto
Date: February 10, 2019
"""
import sys
import argparse
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
sys.path.append("..")
from csa_raster_load import CSASite
# Parse directory name
parser = argparse.ArgumentParser()
parser.add_argument("-d", "--directory",
help="Path to CSA raster data directory")
args = parser.parse_args()
if __name__=="__main__":
site = CSASite(args.directory)
site_width = site.extent_utm[1] - site.extent_utm[0]
site_height = site.extent_utm[3] - site.extent_utm[2]
# Show site info
print("Site width [m]: ", site_width)
print("Site height [m]: ", site_height)
print("Site mosaic res y [m/px]: ", site_height/site.mosaic_rgb.shape[0])
print("Site mosaic res x [m/px]: ", site_width/site.mosaic_rgb.shape[1])
print("Site extent: ", site.extent_utm)
# Plot all maps
fig, axes = plt.subplots(2, 2, figsize=(15,10))
im00 = axes[0,0].imshow(site.mosaic_rgb, extent=site.extent_utm)
axes[0,0].set_title("RGB mosaic")
axes[0,0].set_xlabel('Easting [m]')
axes[0,0].set_ylabel('Northing [m]')
im01 = axes[0,1].imshow(site.dem, extent=site.extent_utm, cmap='gray')
axes[0,1].set_title("Elevation above origin (NW corner) [meters]")
axes[0,1].set_xlabel('Easting [m]')
axes[0,1].set_ylabel('Northing [m]')
im10 = axes[1,0].imshow(site.slope, extent=site.extent_utm, cmap='jet')
axes[1,0].set_title("Slope map [radians]")
axes[1,0].set_xlabel('Easting [m]')
axes[1,0].set_ylabel('Northing [m]')
im11 = axes[1,1].imshow(site.aspect, extent=site.extent_utm, cmap='jet')
axes[1,1].set_title("Aspect map [radians]")
axes[1,1].set_xlabel('Easting [m]')
axes[1,1].set_ylabel('Northing [m]')
plt.suptitle("All georeferenced maps in UTM coordinates (UTM zone {})".format(site.zone_utm))
# Show colorbars
div01 = make_axes_locatable(axes[0,1])
cax01 = div01.append_axes("right", size="5%", pad=0.1)
plt.colorbar(im01, cax=cax01)
div10 = make_axes_locatable(axes[1,0])
cax10 = div10.append_axes("right", size="5%", pad=0.1)
plt.colorbar(im10, cax=cax10)
div11 = make_axes_locatable(axes[1,1])
cax11 = div11.append_axes("right", size="5%", pad=0.1)
plt.colorbar(im11, cax=cax11)
plt.show()