-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmkbathy.py
237 lines (204 loc) · 9.89 KB
/
mkbathy.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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
#!/usr/bin/env python
import sys, os, math
from time import sleep
import numpy as np
from pathlib import Path
import pygmt
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import pymeshlab as ml
from tqdm import trange
# Bathymetry Converter - automation script
###################################################################
PREFIX = 'MontereyBay'
SOURCE = 'bathymetry_source/monterey_13_navd88_2012.nc'
# Flag for single tile generation
ONE_TILE_AT_ORIGIN = False
# Mesh simplification level (0 ~ 3)
SIMPLIFICATION_LEVEL = 0
# Range
STARTLON = -121.825
STARTLAT = 36.790
ENDLON = -121.805
ENDLAT = 36.810
# Resolution
DLON = 0.01
DLAT = 0.01
OVERLON = 0.0005
OVERLAT = 0.0005
###################################################################
TEMPLATE_DIR = '/Bathymetry_Converter/'
if ONE_TILE_AT_ORIGIN:
DLON = ENDLON-STARTLON
DLAT = ENDLAT-STARTLAT
OVERLON = 0.0
OVERLAT = 0.0
print("\n")
print("#--------------------------------------------------------#")
print("#---------------- Bathymetry Converter ----------------#")
print("#--------------------------------------------------------#")
print("# Prefix : " + PREFIX)
print("# Source : " + SOURCE)
print("# Range -- Latitude (N) : " + repr(STARTLAT) + " ~ " + repr(ENDLAT))
print("# \_ longitude(E) : " + repr(STARTLON) + " ~ " + repr(ENDLON))
if ONE_TILE_AT_ORIGIN:
print("# ONE_TILE_AT_ORIGIN flag is ON! Size (Lat, Lon) = (" + repr(DLAT) + ", " + repr(DLON) + ")")
else:
print("# Resolution -- Latitude : " + repr(DLAT) + " with " + repr(OVERLAT) + " overlaps")
print("# \_ longitude : " + repr(DLON) + " with " + repr(OVERLON) + " overlaps")
if SIMPLIFICATION_LEVEL==0:
print("# No simplification")
else:
print("# Simplification level = " + repr(SIMPLIFICATION_LEVEL))
print("#--------------------------------------------------------#\n")
nLon = int(math.ceil((ENDLON-OVERLON-STARTLON)/DLON))
nLat = int(math.ceil((ENDLAT-OVERLAT-STARTLAT)/DLAT))
nTot = nLon * nLat
# Get No Data Value (to skip later)
source = gdal.Open(SOURCE)
band = source.GetRasterBand(1)
band.GetStatistics(True, True)
maxHeight = band.GetMaximum()
NoDataValue = band.GetNoDataValue()
print('Source statistics : ')
print(' Maximum Height : ' + repr(band.GetMaximum()))
print(' Minimum Height : ' + repr(band.GetMinimum()))
print(' NoDataValue represented as : ' + repr(band.GetNoDataValue()))
print(' ')
# Make bathymetry directory (if not exist)
Path('bathymetry').mkdir(parents=True, exist_ok=True) # Template directory
Path(PREFIX).mkdir(parents=True, exist_ok=True) # Final product directory
PATH = "bathymetry/"
# Reproject into WGS (EPSG:4326)
print("Initiating... (It may take longer to initiate for large source)")
os.system("PROJ_NETWORK=ON gdalwarp -t_srs 'EPSG:4326' -of GMT " \
+ SOURCE + " " + PATH + PREFIX + '.grd')
print("Starting conversion..")
longitude = STARTLON
latitude = STARTLAT
# Loop longitude
for lon in trange(nLon):
west = longitude
east = longitude + DLON
west_edge = longitude - OVERLON
east_edge = east + OVERLON
# Loop latitude
for lat in range(nLat):
south = latitude
north = latitude + DLAT
south_edge = latitude - OVERLAT
north_edge = north + OVERLAT
# Create output filename
tileName = "R_" + "{:.3f}".format(west) + '_' + "{:.3f}".format(east) \
+ '_' + "{:.3f}".format(south) + '_' + "{:.3f}".format(north)
# Cut and translate into XYZ point cloud (verbose error only)
xyz = pygmt.grd2xyz(PATH + PREFIX + '.grd', output_type='numpy', \
verbose='e', region=[west_edge, east_edge, south_edge, north_edge])
center_lon = (west_edge + east_edge)/2.0
center_lat = (south_edge + north_edge)/2.0
# Transform project coordniate system from WGS(EPSG:4326) to UTM(EPSG:3857)
source = osr.SpatialReference(); source.ImportFromEPSG(4326)
target = osr.SpatialReference(); target.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(source, target)
for pt in xyz:
if pt[2] == NoDataValue or np.isnan(pt[2]):
pt[2] = -9999
point = ogr.CreateGeometryFromWkt( \
"POINT (" + repr(pt[1]) + " " + repr(pt[0]) + ")")
point.Transform(transform)
if ONE_TILE_AT_ORIGIN:
center = ogr.CreateGeometryFromWkt( \
"POINT (" + repr(center_lat) + " " + repr(center_lon) + ")")
center.Transform(transform)
pt[0] = float(point.ExportToWkt().split(" ")[1].split("(")[1]) \
- float(center.ExportToWkt().split(" ")[1].split("(")[1])
pt[1] = float(point.ExportToWkt().split(" ")[2].split(")")[0]) \
- float(center.ExportToWkt().split(" ")[2].split(")")[0])
# pt[2] = float(pt[2]) - float(maxHeight)
else:
pt[0] = point.ExportToWkt().split(" ")[1].split("(")[1]
pt[1] = point.ExportToWkt().split(" ")[2].split(")")[0]
# Save asc file with 'X Y Z' header for meshlab
filename = PATH + PREFIX + "." + tileName + '.epsg3857'
np.savetxt(filename + '.meshlab.asc', xyz, fmt='%10.5f', delimiter=' ')
os.system("sed -i '1 i\X Y Z' " + filename + '.meshlab.asc')
# Save asc (without header) file after sorting
# https://gdal.org/drivers/raster/xyz.html
# starting with GDAL 3.2.1, cells with same X coordinates must
# be placed on consecutive lines. For a same X coordinate value,
# the columns must be organized by increasing or decreasing Y values.
index = np.lexsort((xyz[:,1],xyz[:,0]))
np.savetxt(filename + '.asc', xyz[index], fmt='%10.5f', delimiter=' ')
# Make texture using colormap
gdal.Translate(filename + '.tif', filename + '.asc', options='-of GTiff -a_srs EPSG:3857')
gdal.DEMProcessing(filename + '.color.tif', filename + '.tif', "color-relief", \
colorFilename='color.txt')
texture = gdal.Translate(filename + '.texture.png', filename + '.color.tif', \
options='-of PNG -ot UInt16 -scale 32.53501 767.4913 0 65535')
# Make mesh (obj) with meshlab
ms = ml.MeshSet()
ms.load_new_mesh(filename + '.meshlab.asc', triangulate=True)
ms.apply_filter('taubin_smooth')
ms.apply_filter('invert_faces_orientation')
ms.apply_filter('parametrization_flat_plane', projectionplane='XY')
if SIMPLIFICATION_LEVEL>=1:
ms.apply_filter('twostep_smooth', normalthr=30.0, stepsmoothnum=5)
if SIMPLIFICATION_LEVEL>=2:
ms.apply_filter('simplification_quadric_edge_collapse_decimation', preserveboundary=True, planarquadric=True)
ms.apply_filter('simplification_quadric_edge_collapse_decimation', preserveboundary=True, planarquadric=True)
if SIMPLIFICATION_LEVEL>=3:
ms.apply_filter('simplification_quadric_edge_collapse_decimation', preserveboundary=True, planarquadric=True, targetperc=0.1)
ms.save_current_mesh(filename + '.obj')
sys.stdout.write("\033[F"+"\r") # supress load_new_mesh output
# Construct model object
TEMPLATE_PATH = TEMPLATE_DIR + 'templates'
MODEL_URI = PREFIX + "." + tileName + '.epsg3857'
if ONE_TILE_AT_ORIGIN:
MODEL_DIR = PREFIX
else:
MODEL_DIR = filename
Path(MODEL_DIR + "/mesh").mkdir(parents=True, exist_ok=True)
os.system("cp " + filename + ".obj " + MODEL_DIR + "/mesh/" )
os.system("cp " + filename + ".obj.mtl " + MODEL_DIR + "/mesh/" )
# create model.config
if ONE_TILE_AT_ORIGIN:
os.system("cat " + TEMPLATE_PATH + "/model.config | sed s#MODEL_NAME#" + PREFIX + "#g > " + MODEL_DIR + "/model.config")
else:
os.system("cat " + TEMPLATE_PATH + "/model.config | sed s#MODEL_NAME#" + PREFIX + "." + tileName + "#g > " + MODEL_DIR + "/model.config")
# create model.sdf
os.system("cat " + TEMPLATE_PATH + "/model.sdf > " + MODEL_DIR +"/model.sdf")
if ONE_TILE_AT_ORIGIN:
os.system("sed -i s#MODEL_NAME#" + PREFIX + "#g " + MODEL_DIR + "/model.sdf")
os.system("sed -i s#MODEL_URI#model://" + PREFIX + "/mesh/" + MODEL_URI + ".obj#g " + MODEL_DIR + "/model.sdf")
else:
os.system("sed -i s#MODEL_NAME#" + PREFIX + "." + tileName + "#g " + MODEL_DIR + "/model.sdf")
os.system("sed -i s#MODEL_URI#model://" + MODEL_URI + "/mesh/" + MODEL_URI + ".obj#g " + MODEL_DIR + "/model.sdf")
# Create material
os.system("cp -r " + TEMPLATE_PATH + "/materials " + MODEL_DIR + "/.")
os.system("cp " + filename + ".texture.png " + MODEL_DIR + "/materials/textures/.")
os.system("sed -i s#TEXTURE_NAME#" + MODEL_URI + ".texture.png#g " + MODEL_DIR + "/materials/scripts/texture.material")
if ONE_TILE_AT_ORIGIN:
os.system("sed -i s#MODEL_NAME#" + PREFIX + "#g " + MODEL_DIR + "/materials/scripts/texture.material")
os.system("sed -i s#TEXTURE_URI#model://" + PREFIX + "#g " + MODEL_DIR + "/model.sdf")
else:
os.system("sed -i s#MODEL_NAME#" + PREFIX + "." + tileName + "#g " + MODEL_DIR + "/materials/scripts/texture.material")
os.system("sed -i s#TEXTURE_URI#model://" + MODEL_URI + "#g " + MODEL_DIR + "/model.sdf")
# Relocate final product
if ONE_TILE_AT_ORIGIN:
os.system("cp -r " + MODEL_DIR + "/* " + PREFIX + "/.")
else:
os.system("cp -r " + MODEL_DIR + " " + PREFIX + "/.")
# update latitude
latitude = north
# update longitude
longitude = east
latitude = STARTLAT
# Remove template directory
os.system("rm -r bathymetry")
# Print generation memo
if ONE_TILE_AT_ORIGIN:
print(" Important memo : ")
print(" A tile is lifted up for " + repr(float(maxHeight)) + " from the original height")
# DONE!
print("\n Done.\n")