GISCourseMaterial/src/map/bm_14.py at master · fmetivier/GISCourseMaterial · GitHub

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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
bm_14.py

(i) Stretch the images according to the histgrams obtained using bm_13
(ii) Merge the bands in one multi-layered raster band
(iii) Reproject and crop the raster
(iv) merge to obtain the final image
(v) plot the map.
"""


import sys
# note that this path is  valid for linux debian
# for windows and mac you need to place a valid path
sys.path.append('/usr/lib/python3/dist-packages/osgeo/utils')


from osgeo import gdal
import gdal_merge as gm
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

# get the GeoTiff driver and register it
driver  =  gdal.GetDriverByName('GTiff')
driver.Register()

# (i)
# Stretch the histograms
# 144030

data = gdal.Open("../../Data/map/raster/China/144030/L5144030_03020060731_B10.TIF")
gdal.Translate(destName = "../../Data/map/raster/China/144030/B10_s.TIF", srcDS = data, scaleParams = [[30,  200,  1,  254]], exponents = [0.7])
data = gdal.Open("../../Data/map/raster/China/144030/L5144030_03020060731_B20.TIF")
gdal.Translate("../../Data/map/raster/China/144030/B20_s.TIF", data,  scaleParams = [[10,  120,  1,  254]], exponents = [0.7])
data = gdal.Open("../../Data/map/raster/China/144030/L5144030_03020060731_B30.TIF")
gdal.Translate("../../Data/map/raster/China/144030/B30_s.TIF", data,  scaleParams = [[8,  140,  1,  254]], exponents = [0.7])


# 145030
data = gdal.Open("../../Data/map/raster/China/145030/L5145030_03020060722_B10.TIF")
gdal.Translate(destName = "../../Data/map/raster/China/145030/B10_s.TIF", srcDS = data, scaleParams = [[30,  200,  1,  254]], exponents = [0.7])
data = gdal.Open("../../Data/map/raster/China/145030/L5145030_03020060722_B20.TIF")
gdal.Translate("../../Data/map/raster/China/145030/B20_s.TIF", data,  scaleParams = [[10,  120,  1,  254]], exponents = [0.7])
data = gdal.Open("../../Data/map/raster/China/145030/L5145030_03020060722_B30.TIF")
gdal.Translate("../../Data/map/raster/China/145030/B30_s.TIF", data,  scaleParams = [[8,  140,  1,  254]], exponents = [0.7])

# (ii)
# Merge the bands
#
gm.main(['', '-o',  '../../Data/map/raster/China/144030/144030_rgb.tif', '-separate', '../../Data/map/raster/China/144030/B30_s.TIF',  '../../Data/map/raster/China/144030/B20_s.TIF',  '../../Data/map/raster/China/144030/B10_s.TIF'])
gm.main(['', '-o',  '../../Data/map/raster/China/145030/145030_rgb.tif', '-separate', '../../Data/map/raster/China/145030/B30_s.TIF',  '../../Data/map/raster/China/145030/B20_s.TIF',  '../../Data/map/raster/China/145030/B10_s.TIF'])


# (iii)
# Crop and reproject the images
#
gdal.Warp("../../Data/map/raster/China/144030/144030_rgb_ll.tif", "../../Data/map/raster/China/144030/144030_rgb.tif",  dstSRS = "+proj=latlong +datum=WGS84", outputBounds = (84.25,  42.5,  86.1,  43.5)) # NB: no spaces surrounding '=' in the SRS !!!
gdal.Warp("../../Data/map/raster/China/145030/145030_rgb_ll.tif", "../../Data/map/raster/China/145030/145030_rgb.tif",  dstSRS = "+proj=latlong +datum=WGS84", outputBounds = (82.5,  42.5,  84.5,  43.5)) # NB: no spaces surrounding '=' in the SRS !!!

# (iv)
# Merge into the final image
#
gm.main(['', '-n', '0', '-o',  '../../Data/map/raster/China/bayan.tif',  '../../Data/map/raster/China/144030/144030_rgb_ll.tif',  '../../Data/map/raster/China/145030/145030_rgb_ll.tif'])


# (v)
# Draw the map
#
ds  =  gdal.Open("../../Data/map/raster/China/bayan.tif")
data  =  ds.ReadAsArray()


fig = plt.figure(figsize = (20, 10))# because it is a raster the output size is important
map = Basemap(projection = 'cyl',  llcrnrlon = 82.5,  llcrnrlat = 42.5, urcrnrlon = 86.1, urcrnrlat = 43.5)

map.imshow(plt.imread("../../Data/map/raster/China/bayan.tif"),  origin = 'upper')

map.drawmeridians(range(83, 88, 1), color = 'r', labels = [0, 0, 1, 1], fontsize = 14)
map.drawparallels([42.5, 43, 43.5], color = 'r', labels = [1, 1, 0, 0], fontsize = 14)

plt.show()