/
ag_build_elevation_map.py
78 lines (52 loc) · 2.44 KB
/
ag_build_elevation_map.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
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 17 18:24:29 2019
@author: Ethan
"""
def buildElevationMap(latPixels=360, lonPixels=720):
import numpy as np
import rasterio
import glob
import sys
import pickle
dataDir = '/home/edcoffel/drive/MAX-Filer/Research/Climate-01/Personal-F20/edcoffel-F20/data/elevation'
files = glob.glob('%s/*.tif'%dataDir)
lat = np.linspace(90, -90, latPixels)
lon = np.linspace(-180, 180, lonPixels)
lonStep = lon[1]-lon[0]
latStep = lat[1]-lat[0]
elevationMap = np.zeros([len(lat), len(lon)])
for file in files:
print('processing %s...'%file)
elev = rasterio.open(file)
elevData = elev.read(1)
tl_lat = elev.bounds.top
tl_lon = elev.bounds.left
br_lat = elev.bounds.bottom
br_lon = elev.bounds.right
xStep = elev.transform[0]
yStep = elev.transform[4]
mapYLeft = np.where(abs(lon-tl_lon) == np.nanmin(abs(lon-tl_lon)))[0][0]
mapYRight = np.where(abs(lon-br_lon) == np.nanmin(abs(lon-br_lon)))[0][0]
mapXTop = np.where(abs(lat-tl_lat) == np.nanmin(abs(lat-tl_lat)))[0][0]
mapXBottom = np.where(abs(lat-br_lat) == np.nanmin(abs(lat-br_lat)))[0][0]
# print((mapYLeft, mapYRight, mapXTop, mapXBottom))
for x in range(min(mapXTop, mapXBottom), max(mapXTop, mapXBottom)+1):
for y in range(min(mapYLeft, mapYRight), max(mapYLeft, mapYRight)+1):
if x >= elevationMap.shape[0] or y >= elevationMap.shape[1]:
continue
mapCoordX_tl = int(round((tl_lat-lat[x])/xStep))
mapCoordX_br = int(round((tl_lat-(lat[x]+latStep))/xStep))
mapCoordY_tl = int(round((tl_lon-lon[y])/yStep))
mapCoordY_br = int(round((tl_lon-(lon[y]+lonStep))/yStep))
if mapCoordY_br >= elevData.shape[1]:
mapCoordY_br = elevData.shape[1]
if mapCoordX_br >= elevData.shape[0]:
mapCoordX_br = elevData.shape[0]
meanElev = np.nanmean(np.nanmean(elevData[mapCoordX_tl:mapCoordX_br, mapCoordY_tl:mapCoordY_br]))
if not np.isnan(meanElev):
elevationMap[x,y] = meanElev
with open('%s/elevation-map-%d-%d.dat'%(dataDir, latPixels, lonPixels), 'wb') as f:
e = {'lat':lat, 'lon':lon, 'elevationMap':elevationMap}
pickle.dump(e, f)
return lat, lon, elevationMap