hycom_uv_netcdf_to_veclocity_gtiff.py

#!/usr/bin/env python

# Warning: This is not great code.

# Copyright 2015 Google Inc. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http:#www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Draft script to convert HYCOM netcdf to gtiff.
#
# Convert the two component speed files to velocity.
# Need to add direction using hypot().
#
# This script needs serious cleanup.
#
# Requires that you start with the regular grid versions of HYCOM uv velocity
# files.  e.g.
#
#   ftp://ftp.hycom.org/datasets/GLBu0.08/expt_91.1/data/hindcasts/2015/
#     hycom_glb_911_2015093000_t000_uv3z.nc
#
# Author: Kurt Schwehr

import math
import sys

import numpy
from osgeo import gdal
from osgeo import osr

OLD_NODATA = -30000
NEW_NODATA = -1
OLD_SCALE = 0.001

def Dist(a, b):
  if a == OLD_NODATA or b == OLD_NODATA:
    return NEW_NODATA  # Redefine
  #  Also add the scale factor

  # Convert to actual meters per second.
  return OLD_SCALE * math.sqrt(a*a + b*b)

VectorizedDist = numpy.vectorize(Dist)

def ConvertLayers(src_filename, dst_filename):
  src_filename_u = 'NETCDF:%s:water_u' % src_filename
  src_filename_v = 'NETCDF:%s:water_v' % src_filename

  src_u = gdal.Open(src_filename_u, gdal.GA_ReadOnly)
  src_v = gdal.Open(src_filename_v, gdal.GA_ReadOnly)

  w = src_u.RasterXSize
  h = src_u.RasterYSize

  num_bands = src_u.RasterCount
  if num_bands != 40:
    printf('Error: num_bands != 40')
    sys.exit(1)

  # dst = gdal.Open(dst_filename, gdal.GA_Update)
  tiff_drv = gdal.GetDriverByName('GTiff')

  # options = ['compress=deflate', 'predictor=2']
  options = []
  data_type = gdal.GDT_Float32

  dst = tiff_drv.Create(dst_filename, w, h, num_bands, data_type, options)

  gt = src_u.GetGeoTransform()
  dst_gt = (gt[0] - 180, gt[1], gt[2], gt[3], gt[4], gt[5])
  dst.SetGeoTransform(dst_gt)

  dst.SetProjection(osr.SRS_WKT_WGS84)

  metadata = src_u.GetMetadata()
  metadata.update(src_v.GetMetadata())
  dst.SetMetadata(metadata)

  for band_num in range(1, num_bands):
    print 'band num:', band_num

    band_u = src_u.GetRasterBand(band_num)
    band_v = src_v.GetRasterBand(band_num)

    band_metadata = band_u.GetMetadata()
    band_metadata['NETCDF_VARNAME'] = 'speed'
    band_metadata['_FillValue'] = str(NEW_NODATA)
    band_metadata['long_name'] = 'Sea Water Speed'
    band_metadata['standard_name'] = 'sea_water_speed'
    band_metadata['scale_factor'] = '1'

    band_dst = dst.GetRasterBand(band_num)
    band_dst.SetMetadata(band_metadata)
    band_dst.SetNoDataValue(NEW_NODATA)
    for row in range(h):
      if row % 200 == 0:
        print 'row', row
      line_u = band_u.ReadAsArray(0, row, w, 1, w, 1)[0].astype(numpy.float)
      line_v = band_v.ReadAsArray(0, row, w, 1, w, 1)[0].astype(numpy.float)
      speed = VectorizedDist(line_u, line_v)
      speed = numpy.concatenate(numpy.hsplit(speed, 2)[::-1])
      band_dst.WriteArray(speed.reshape((1, w)), 0, row)

  tiff_drv = None

def main():
  ConvertLayers('hycom_glb_911_2015102200_t000_uv3z.nc',
                'hycom_glb_911_2015102200_t000_uv3z-speed.tif')

if __name__ == '__main__':
    main()