/
hillshade_mem.py
152 lines (131 loc) · 4.93 KB
/
hillshade_mem.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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst
import sys
import os
import pdb
## Crop a global epsg:4326 DEM in 10° squares, compute hillshading and
## wrap output to epsg:3857. No edge effects
## takes ~9h on 8 cores, nvme disk for a 1 arcsec / 30m DEM
gdal.SetConfigOption('GDAL_NUM_THREADS', '8') # speedup warping
infile=sys.argv[1]
input_DEM=gdal.Open( infile, gdal.GA_ReadOnly )
src_res=0.0002777777777778966603
# ~ for lat in range(80,90,10):
# ~ for lon in range(-180,180,10):
# ~ for lat in range(-70,70,10):
# ~ for lon in range(-180,180,10):
# ~ ulx = lon-2*src_res
# ~ uly = lat+10+2*src_res
# ~ lrx= lon+10+2*src_res
# ~ lry= lat-2*src_res
# ~ if ulx <= -180: ulx =-179.999
# ~ if lrx >= 180: lrx=179.999
# ~ if uly >= 90: uly=89.999
# ~ if lry<=-90: lry=-89.999
# ~ fulx = lon
# ~ fuly = lat+10
# ~ flrx= lon+10
# ~ flry= lat
# ~ if fulx <= -180: fulx =-179.999
# ~ if flrx >= 180: flrx=179.999
# ~ if fuly >= 90: fuly=89.999
# ~ if flry <=-90: flry=-89.999
# ~ proj_win=[ulx,uly,lrx,lry]
# ~ final_win=[fulx,fuly,flrx,flry]
# ~ cropped_DEM = gdal.GetDriverByName("Memory").Create("cropped_DEM", 0, 0, 0, gdal.GDT_Unknown)
# ~ tmp_Hillshade = gdal.GetDriverByName("Memory").Create("tmp_Hillshade", 0, 0, 0, gdal.GDT_Unknown)
# ~ tmp_Hillshade_3857 = gdal.GetDriverByName("Memory").Create("tmp_Hillshade_3857", 0, 0, 0, gdal.GDT_Unknown)
# ~ final_Hillshade="../hillshade_tiles/mem_HS_tile_3857_"+str(lon)+"_"+str(lat)+".tif"
# ~ print(proj_win, final_win)
# ~ print("crop...")
# ~ gdal.Translate("/vsimem/cropped_DEM", input_DEM, projWin=proj_win)
# ~ print("hillshading...")
# ~ gdal.DEMProcessing("/vsimem/tmp_Hillshade",
# ~ "/vsimem/cropped_DEM",
# ~ "hillshade",
# ~ zFactor= 2,
# ~ scale=111120,
# ~ azimuth= 320.0,
# ~ igor=True)
# ~ print("warping...")
# ~ gdal.Warp("/vsimem/tmp_Hillshade_3857",
# ~ "/vsimem/tmp_Hillshade",
# ~ srcSRS="epsg:4326",
# ~ dstSRS="epsg:3857",
# ~ xRes=30,
# ~ yRes=30,
# ~ resampleAlg="cubic")
# ~ print("shaving...")
# ~ gdal.Translate(final_Hillshade,"/vsimem/tmp_Hillshade_3857",
# ~ projWin=final_win,
# ~ projWinSRS="epsg:4326",
# ~ format="GTiff",
# ~ creationOptions=["TILED=YES","NUM_THREADS=4","BIGTIFF=YES","COMPRESS=DEFLATE"])
# creationOptions=["TILED=YES","NUM_THREADS=4","BIGTIFF=YES","COMPRESS=ZSTD","ZSTD_LEVEL=2"]) ZSTD may cause incompativbility with mapnik
# ~ try:
# ~ os.remove(cropped_DEM)
# ~ os.remove(tmp_Hillshade)
# ~ os.remove(tmp_Hillshade_3857)
# ~ except: pass
# ~ print("done.\n")
for lat in range(70,72,2):
for lon in range(-180,180,10):
ulx = lon-2*src_res
uly = lat+2+2*src_res
lrx= lon+10+2*src_res
lry= lat-2*src_res
if ulx <= -180: ulx =-179.999
if lrx >= 180: lrx=179.999
if uly >= 90: uly=89.999
if lry<=-90: lry=-89.999
fulx = lon
fuly = lat+2
flrx= lon+10
flry= lat
if fulx <= -180: fulx =-179.999
if flrx >= 180: flrx=179.999
if fuly >= 90: fuly=89.999
if flry <=-90: flry=-89.999
proj_win=[ulx,uly,lrx,lry]
final_win=[fulx,fuly,flrx,flry]
cropped_DEM = gdal.GetDriverByName("Memory").Create("cropped_DEM", 0, 0, 0, gdal.GDT_Unknown)
tmp_Hillshade = gdal.GetDriverByName("Memory").Create("tmp_Hillshade", 0, 0, 0, gdal.GDT_Unknown)
tmp_Hillshade_3857 = gdal.GetDriverByName("Memory").Create("tmp_Hillshade_3857", 0, 0, 0, gdal.GDT_Unknown)
final_Hillshade="../hillshade_tiles/mem_HS_tile_3857_"+str(lon)+"_"+str(lat)+".tif"
print(proj_win, final_win)
print("crop...")
gdal.Translate("/vsimem/cropped_DEM", input_DEM, projWin=proj_win)
print("hillshading...")
gdal.DEMProcessing("/vsimem/tmp_Hillshade",
"/vsimem/cropped_DEM",
"hillshade",
zFactor= 2,
scale=111120,
azimuth= 320.0,
igor=True)
print("warping...")
gdal.Warp("/vsimem/tmp_Hillshade_3857",
"/vsimem/tmp_Hillshade",
srcSRS="epsg:4326",
dstSRS="epsg:3857",
xRes=30,
yRes=30,
resampleAlg="cubic")
print("shaving...")
gdal.Translate(final_Hillshade,"/vsimem/tmp_Hillshade_3857",
projWin=final_win,
projWinSRS="epsg:4326",
format="GTiff",
creationOptions=["TILED=YES","NUM_THREADS=4","BIGTIFF=YES","COMPRESS=DEFLATE"])
# ~ creationOptions=["TILED=YES","NUM_THREADS=4","BIGTIFF=YES","COMPRESS=ZSTD","ZSTD_LEVEL=2"]) ZSTD may cause incompativbility with mapnik
try:
os.remove(cropped_DEM)
os.remove(tmp_Hillshade)
os.remove(tmp_Hillshade_3857)
except: pass
print("done.\n")
# ~ gdal_translate --config GDAL_CACHEMAX 1000 --config GDAL_DISABLE_READDIR_ON_OPEN TRUE -co "BIGTIFF=YES" -co "TILED=YES" -co "COMPRESS=DEFLATE" -co "NUM_THREADS=8" ../hillshade.vrt ../Hillshade3.tif
# ~ gdaladdo -ro --config COMPRESS_OVERVIEW DEFLATE --config GDAL_NUM_THREADS 8 --config GDAL_CACHEMAX 1000 Hillshade3.tif 2 4 8 16 32