Skip to content

Commit

Permalink
Merge pull request #836 from mnagaso/devel
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielpeter
danielpeter committed Apr 15, 2024
2 parents fefc7ee + f5d709e commit 1642be7
Showing 1 changed file with 34 additions and 15 deletions.
49 changes: 34 additions & 15 deletions DATA/topo_bathy/run_create_topo_bathy_file.py
View file
Expand Up @@ -53,6 +53,19 @@
if gmt_major >= 5 and gmt_minor >= 3: print("version ok")
print("")

# xyz2grd, grd2xyz, grdsample, grdinfo is wrapped in gmt after version 6.0
if gmt_major >= 6:
xyz2grd = 'gmt xyz2grd'
grd2xyz = 'gmt grd2xyz'
grdsample = 'gmt grdsample'
grdinfo = 'gmt grdinfo'
gmtinfo = 'gmt info'
else:
xyz2grd = 'xyz2grd'
grd2xyz = 'grd2xyz'
grdsample = 'grdsample'
grdinfo = 'grdinfo'
gmtinfo = 'gmtinfo'

####################################################################
# User parameters
Expand Down Expand Up @@ -362,21 +375,21 @@ def resample_topo_file(topo,filename_web,filename_grid):

if topo == 'etopo1':
# converting from raw binary to grid file
cmd = 'xyz2grd ' + filename_web + ' -V -Rd ' + grid + ' -r -ZTLh ' + interval_min
cmd = f'{xyz2grd} ' + filename_web + ' -V -Rd ' + grid + ' -r -ZTLh ' + interval_min

elif topo == 'etopo2' or topo == 'etopo4' or topo == 'etopo5' or topo == 'etopo15':
if filename_web == 'ETOPO2v2c_i2_LSB.bin':
# converting from etopo2v2c raw binary to 2m grid file first
print(" converting binary file to grid file: etopo2v2c.grd")
cmd = 'xyz2grd ETOPO2v2c_i2_LSB.bin -Rd -I2m -Getopo2v2c.grd -r -ZTLh -V'
cmd = f'{xyz2grd} ETOPO2v2c_i2_LSB.bin -Rd -I2m -Getopo2v2c.grd -r -ZTLh -V'
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
# resample
cmd = 'grdsample etopo2v2c.grd -V -Rd ' + grid + ' ' + interval_min
cmd = f'{grdsample} etopo2v2c.grd -V -Rd ' + grid + ' ' + interval_min
else:
# direct resample
cmd = 'grdsample ' + filename_web + ' -V -Rd ' + grid + ' ' + interval_min
cmd = f'{grdsample} ' + filename_web + ' -V -Rd ' + grid + ' ' + interval_min

elif 'mars' in topo:
# converting to raw binary grid file
Expand All @@ -391,12 +404,12 @@ def resample_topo_file(topo,filename_web,filename_grid):
# > grdimage topo.grd -JX6i+ -I+d -P -Vl > topo.ps
print(" converting binary file to grid file: "+filename_web)
# using int-16 values
cmd = 'xyz2grd ' + filename_web + ' -Gmarstopo1.875.grd -R0/360/-90/90 -I1.875m -r -ZTLhw -fic -fog -V'
cmd = f'{xyz2grd} ' + filename_web + ' -Gmarstopo1.875.grd -R0/360/-90/90 -I1.875m -r -ZTLhw -fic -fog -V'
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
# resample
cmd = 'grdsample marstopo1.875.grd -V -Rd ' + grid + ' ' + interval_min
cmd = f'{grdsample} marstopo1.875.grd -V -Rd ' + grid + ' ' + interval_min
elif 'moon' in topo:
# converting to raw binary grid file
# 0.1 degree resolution = 0.1 / (1.60.) = 6 arc-minute
Expand Down Expand Up @@ -441,27 +454,27 @@ def resample_topo_file(topo,filename_web,filename_grid):
# workaround: directly using xyz2grd fails, complains about wrong array sizes, probably due to some header data.
# thus, we output just the elevation data and re-read into grd-format
# converts to binary float without x,y positions
cmd = 'grd2xyz topo.grd -ZTLf > topo.bin'
cmd = f'{grd2xyz} topo.grd -ZTLf > topo.bin'
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
# using int-16 values from 16-pixel resolution = 0.0625 degree resolution
if topo == 'moon4':
cmd = 'xyz2grd topo.bin -Gtopo.moon.grd -R0/360/-90/90 -I0.0625 -r -ZTLf -fic -fog -Vl'
cmd = f'{xyz2grd} topo.bin -Gtopo.moon.grd -R0/360/-90/90 -I0.0625 -r -ZTLf -fic -fog -Vl'
else:
# using int-16 values from 64-pixel resolution = 0.9375m arc-minutes
cmd = 'xyz2grd topo.bin -Gtopo.moon.grd -R0/360/-90/90 -I0.9375m -r -ZTLf -fic -fog -Vl'
cmd = f'{xyz2grd} topo.bin -Gtopo.moon.grd -R0/360/-90/90 -I0.9375m -r -ZTLf -fic -fog -Vl'
#cmd = 'xyz2grd ' + filename_web + ' -Gmoontopo.grd -R0/360/-90/90 -I0.9375m -r -ZTLhw -fic -fog -V'
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
print("")
cmd = 'grdinfo topo.moon.grd'
cmd = f'{grdinfo} topo.moon.grd'
status = subprocess.call(cmd, shell=True)
check_status(status)
print("")
# resample
cmd = 'grdsample topo.moon.grd -V -Rd ' + grid + ' ' + interval_min
cmd = f'{grdsample} topo.moon.grd -V -Rd ' + grid + ' ' + interval_min

else:
print("invalid topo %s for resampling" % topo)
Expand All @@ -486,23 +499,23 @@ def extract_to_ascii(topo,filename_in,filename_out):
print(" file %s exists already, skipping extraction" % filename_out)
else:
# file info
cmd = 'grdinfo ' + filename_in
cmd = f'{grdinfo} ' + filename_in
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
print("")

# converts to ascii, extracts elevation values
print(" extracting elevation values to ascii file...")
cmd = 'grd2xyz -Rg ' + filename_in + ' | awk \'{ print $3 }\' > ' + filename_out
cmd = f'{grd2xyz} -Rg ' + filename_in + ' | awk \'{ print $3 }\' > ' + filename_out
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
print("written to: %s" % filename_out)

print("")
print(" gmt info:")
cmd = 'gmtinfo ' + filename_out
cmd = f'{gmtinfo} ' + filename_out
print("> ",cmd)
status = subprocess.call(cmd, shell=True)
check_status(status)
Expand Down Expand Up @@ -640,7 +653,13 @@ def smooth_topo_bathy(topo,filename_in,filename_out):
# sys.exit(1)
print(" size of window filter is (%d x %d)" %(2*SIZE_FILTER_ONE_SIDE+1,2*SIZE_FILTER_ONE_SIDE+1))
print("")
area_window = np.float((2*SIZE_FILTER_ONE_SIDE+1)**2)
area_window = (2*SIZE_FILTER_ONE_SIDE+1)**2

# Check numpy version
if np.__version__ >= '1.24':
area_window = float(area_window)
else:
area_window = np.float(area_window)


count = 0
Expand Down

0 comments on commit 1642be7

Please sign in to comment.