hycom_uv_netcdf_to_uv_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)

def Scale(a):
  if a == OLD_NODATA:
    return NEW_NODATA  # Redefine

  # Convert to actual meters per second.
  return OLD_SCALE * a

VectorizedDist = numpy.vectorize(Dist)
VectorizedScale = numpy.vectorize(Scale)

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])
  # print 'GeoTransform:', dst_gt
  dst.SetGeoTransform(dst_gt)

  dst.SetProjection(osr.SRS_WKT_WGS84)

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

  out_band = 0
  for band_num in range(1, 21):
    print 'band num:', band_num

    band_u = src_u.GetRasterBand(band_num)

    band_metadata = band_u.GetMetadata()
    band_metadata['NETCDF_VARNAME'] = 'water_u'
    band_metadata['_FillValue'] = str(NEW_NODATA)
    band_metadata['long_name'] = 'Eastward Water Velocity'
    band_metadata['standard_name'] = 'eastward_sea_water_velocity'
    band_metadata['scale_factor'] = '1'

    out_band = out_band + 1
    band_dst = dst.GetRasterBand(out_band)
    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)
      udata = VectorizedScale(line_u)
      udata = numpy.concatenate(numpy.hsplit(udata, 2)[::-1])
      band_dst.WriteArray(udata.reshape((1, w)), 0, row)

    band_v = src_v.GetRasterBand(band_num)

    band_metadata = band_v.GetMetadata()
    band_metadata['NETCDF_VARNAME'] = 'water_v'
    band_metadata['_FillValue'] = str(NEW_NODATA)
    band_metadata['long_name'] = 'Northward Water Velocity'
    band_metadata['standard_name'] = 'northward_sea_water_velocity'
    band_metadata['scale_factor'] = '1'

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

  tiff_drv = None

def main():
  ConvertLayers('hycom_glb_912_2016060500_t000_uv3z.nc',
                'hycom_glb_912_2016060500_t000_uv3z.tif')

if __name__ == '__main__':
    main()

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()

mapfile-for-hycom-ocean-currents.map

MAP

# ... add appropriate header info here

SYMBOL
  NAME "horizline"
  TYPE VECTOR
  POINTS
     0 0
     1 0
  END # points
END # symbol
SYMBOL
  NAME "arrowhead"
  TYPE vector
  FILLED true
  #ANCHORPOINT 0 0.5
  POINTS
    0 2
    4 1
    0 0
  END # points
END # symbol
SYMBOL
  NAME "arrowtail"
  TYPE vector
  FILLED true
  ANCHORPOINT 1 0.5 # to shift the arrowtail
  POINTS
    0 2
    4 1
    0 0
    -99 -99
    0 1
    4 1
  END # points
END # symbol
LAYER
  NAME "my_uv_test"
  TYPE POINT
  STATUS DEFAULT
  CONNECTIONTYPE uvraster
  DATA "hycom_glb_912_2016060500_t000_uv3z.tif"
  PROCESSING "BANDS=1,2"
#  PROCESSING "BANDS=9,10"
  PROCESSING "UV_SPACING=20"
  PROCESSING "UV_SIZE_SCALE=30.0"
  CLASS
    STYLE
      SYMBOL "horizline"
      ANGLE [uv_angle]
      SIZE [uv_length]
      WIDTH 3
      #COLOR 100 255 0
      COLORRANGE 0 255 0 255 0 0
      DATARANGE 0 25
      RANGEITEM "uv_length"
    END # style
    STYLE
      SYMBOL "arrowhead"
      ANGLE [uv_angle]
      SIZE 7
      #COLOR 255 0 0
      COLORRANGE 0 255 0 255 0 0
      DATARANGE 0 25
      RANGEITEM "uv_length"
      POLAROFFSET [uv_length_2] [uv_angle]
  END # class
END # layer

END # map