Skip to content

Instantly share code, notes, and snippets.

@grinsted
Created June 23, 2020 07:51
Show Gist options
  • Save grinsted/b5d365853f98525f730722c7e1035a65 to your computer and use it in GitHub Desktop.
Save grinsted/b5d365853f98525f730722c7e1035a65 to your computer and use it in GitHub Desktop.
import subprocess
import os
import glob
for filename in glob.iglob(os.path.join(".", "**", "*.nc"), recursive=True):
if os.stat(filename).st_size < 50:
continue
result = subprocess.run(["gdalinfo.exe", filename], stdout=subprocess.PIPE)
if result.stdout.rfind(b"COMPRESSION=") > 0:
continue
print("\n\n")
print(filename)
directions = {
"vx": "land_ice_surface_easting_velocity",
"vy": "land_ice_surface_northing_velocity",
}
for key, value in directions.items():
vrtfilename = filename.replace(".nc", "_{}.vrt".format(key))
outputfilename = filename.replace(".nc", "_{}.tif".format(key))
# MAKE A VRT file using the netcdf driver:
command = [
"gdal_translate.exe",
"-a_srs",
"EPSG:3413",
"NETCDF:{}:{}".format(filename, value),
vrtfilename,
]
result = subprocess.run(command, stdout=subprocess.PIPE)
# MAGIC replace in vrt file:
with open(vrtfilename, "r") as file:
vrtfiledata = file.read()
vrtfiledata = vrtfiledata.replace("NETCDF:", "HDF5:")
vrtfiledata = vrtfiledata.replace(":land_ice_surface", "://land_ice_surface")
# Write the file out again
with open(vrtfilename, "w") as file:
file.write(vrtfiledata)
# do the conversion to TIF using the HDF5 driver
command = [
"gdal_translate.exe",
"-a_srs",
"EPSG:3413",
"-co",
"COMPRESS=DEFLATE",
"-co",
"predictor=2",
"-co",
"tiled=yes",
vrtfilename,
outputfilename,
]
result = subprocess.run(command, stdout=subprocess.PIPE)
print(result.stdout.decode())
if result.stdout.rfind(b"done") < 0:
print("WHAT?")
break
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment