Fancy NetCDF Python Example

download_hycom.py

# Just an example it wont outside of the application is was built in
# Uses Scipy.io netcdf and 'snetCDF 4 Dataset
from scipy.io import netcdf
from netCDF4 import Dataset


def navy_hycom_download(level='top'):
      from pl_plot.plotter import HycomPlotter
      print "NAVY HYCOM DOWNLOAD"

      # TODO: Interpolate

      top = 0.0
      bot = 2500.0
      vert_coord= str(top)

      time_start_dt = DataFileManager.get_last_forecast_for_roms()
      time_start = create_nomads_string_from_datetime(DataFileManager.get_last_forecast_for_roms())
      time_end = create_nomads_string_from_datetime(DataFileManager.get_last_forecast_for_roms() + timedelta(days=30)) # Essentially set to infinity


      URL = 'http://ncss.hycom.org/thredds/ncss/GLBv0.08/expt_93.0/data/forecasts/Forecast_Mode_Run_(FMRC)_best.ncd?' \
            'var=salinity&var=water_temp&var=water_u&var=water_v' \
            '&north=50&west=-150&east=-110&south=25&disableProjSubset=on' \
            '&horizStride=1' \
            '&time_start='+time_start+'' \
            '&time_end='+time_end+'' \
            '&timeStride=1' \
            '&vertCoord='+vert_coord+'' \
            '&addLatLon=true&accept=netcdf'

      print "DOWNLOADING HYCOM DF...",
      destination_directory = os.path.join(settings.MEDIA_ROOT, settings.NETCDF_STORAGE_DIR)
      local_filename = "{0}_{1}.nc".format(settings.NAVY_HYCOM_DF_FN, create_string_from_dt_for_file(time_start_dt))
      urllib.urlretrieve(url=URL, filename=os.path.join(destination_directory, local_filename))
      print "FILE DOWNLOADED...",

      # Application specific wrapper to open the file we just downloaded. It just calls
      # netcdf.netcdf_file(..)
      plotter = HycomPlotter(local_filename)

      time='time'
      salt='salinity'
      temp='water_temp
      cur_u='water_u'
      cur_v='water_v'

      depth = plotter.data_file.variables['depth']
      times = plotter.data_file.variables[time]
      salinity = plotter.data_file.variables[salt][:, 0, :, :]
      temps = plotter.data_file.variables[temp][:, 0, :, :]
      curs_u = plotter.data_file.variables[cur_u][:, 0, :, :]
      curs_v = plotter.data_file.variables[cur_v][:, 0, :, :]

      # Squeeze out the depth of each variable...
      salinity = numpy.squeeze(salinity)
      temps = numpy.squeeze(temps)
      curs_u = numpy.squeeze(curs_u)
      curs_v = numpy.squeeze(curs_v)

      lons = plotter.data_file.variables['lon']
      lats = plotter.data_file.variables['lat']

      basetime = datetime.strptime(plotter.data_file.variables[time].units, 'hours since %Y-%m-%d 12:00:00.000 UTC')

      ts1 = times; ts1 = map(float, ts1)

      dates = [(basetime + n * timedelta(hours=4)) for n in range(times.shape[0])]

      ts2 = date2num(dates[:], units=times.units, calendar=times.calendar)

      for i in range(len(ts2)):
          ts2[i] += times[0] + 1 # Now bump up the times to be equal with what we desire


      salinity_int = numpy.empty([ts2.shape[0], lats.shape[0], lons.shape[0]]) # Array to be filled
      temps_int = numpy.empty([ts2.shape[0], lats.shape[0], lons.shape[0]]) # Array to be filled
      curs_u_int = numpy.empty([ts2.shape[0], lats.shape[0], lons.shape[0]]) # Array to be filled
      curs_v_int = numpy.empty([ts2.shape[0], lats.shape[0], lons.shape[0]]) # Array to be filled

      print "...Interpolating...",

      for i in range(lats.shape[0]):
          for j in range(lons.shape[0]):
              """ For each lat and lon across t, time, interpolate from ts1, orginal timestamp,
              the the new timestamp, ts2.
              """
              salinity_int[:, i, j] = numpy.interp(ts2, ts1, salinity[:, i, j]) # Heights interpolated
              temps_int[:, i, j] = numpy.interp(ts2, ts1, temps[:, i, j])
              curs_u_int[:, i, j] = numpy.interp(ts2, ts1, curs_u[:, i, j])
              curs_v_int[:, i, j] = numpy.interp(ts2, ts1, curs_v[:, i, j])

      plotter.close_file() # Close readonly Datafile
      print "saving interpolated values..."

      # Another wrapper to open the netcdf file.
      # Uses the Dataset call from netcdf. See above
      plotter.write_file(local_filename) # Open datafile for writing

      # Write values
      plotter.data_file.variables[time][:] = ts2
      plotter.data_file.variables[salt][:, :, :] = salinity_int
      plotter.data_file.variables[temp][:, :, :] = temps_int
      plotter.data_file.variables[cur_u][:, :, :] = curs_u_int
      plotter.data_file.variables[cur_v][:, :, :] = curs_v_int

      plotter.close_file() # Close rw datafile

      datafile = DataFile(
          type='HYCOM',
          download_datetime=timezone.now(),
          generated_datetime=timezone.now(),
          model_date=time_start_dt,
          file=local_filename,
      )
      datafile.save() # S-S-Save that datafile entry!
      print "finished! Datafile saved!"

return datafile.id