GISCourseMaterial/src/map/cm_12_2.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
84
85
86
87
88
89
90
91
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
cm_12-2.py


Plot a georeferenced DEM
using merging and cropping

Cartopy version
"""

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from osgeo import gdal
from osgeo import osr
import pyproj
from numpy import linspace
from numpy import meshgrid
from numpy import floor, ceil


fname = ["../../Data/map/raster/Guadeloupe/ASTGTM2_N16W062/ASTGTM2_N16W062_dem.tif",\
"../../Data/map/raster/Guadeloupe/ASTGTM2_N15W062/ASTGTM2_N15W062_dem.tif"]

fig = plt.figure()
ax = fig.add_subplot(111, projection = ccrs.PlateCarree())

minx=180
maxx=-180
miny=90
maxy=-90

for f in fname:

	ds  =  gdal.Open(f)
	data  =  ds.ReadAsArray()


	# map limits
	geotransform  =  ds.GetGeoTransform()
	proj=ds.GetProjection()
	print(proj)
	originX  =  geotransform[0]
	originY  =  geotransform[3]
	print( originX,  originY )
	pixelWidth  =  geotransform[1]
	pixelHeight  =  geotransform[5]
	print( pixelWidth,  pixelHeight )

	w = abs(pixelWidth * ds.RasterXSize)
	h = abs(pixelHeight * ds.RasterYSize)

	#corners
	lx = originX
	ly = originY-h
	rx = lx+w
	ry = originY

	if minx > min(lx,rx):
		minx = min(lx,rx)
	if maxx < max(lx,rx):
		maxx = max(lx,rx)
	if miny > min(ly,ry):
		miny = min(ly,ry)
	if maxy < max(ly,ry):
		maxy = max(ry,ly)

	extent = (originX, originX + abs(pixelWidth * ds.RasterXSize)  , originY - abs(pixelHeight * ds.RasterYSize) , originY)
	print(extent)

	print('lower left corner', lx, ly )
	print('upper right corner',  rx, ry )

	#get the spatial projection of the DEM

	inproj = osr.SpatialReference()
	inproj.ImportFromWkt(proj)
	projcs = inproj.GetAuthorityCode('PROJCS')
	print(projcs)
	#~ DEM_proj = ccrs.epsg(6326)


	ax.imshow(data, extent=extent, origin='upper', cmap='terrain', transform=ccrs.PlateCarree())


ax.coastlines(resolution='10m')
ax.set_xlim(minx,maxx)
ax.set_ylim(miny,maxy)
# plt.savefig("../../figures/antilles.svg", bbox_inches = "tight")
plt.show()