module AgwxBiophys::DegreeDays

Constants

DD_MAX_TEMP
M_1_PI

Public Instance Methods

build_averages(min_grids,max_grids,start_year=START_YEAR,end_year=END_YEAR) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 149
def build_averages(min_grids,max_grids,start_year=START_YEAR,end_year=END_YEAR)
  averages = {}
  layers_for_doy = {}
  # print "building averages. "
  # check how many years have data for each doy
  (START_DOY..END_DOY).each do |doy|
    for year in (start_year..end_year)
      # Does this year have this DOY?
      if min_grids[year].get_by_index(0,0,doy)
        if layers_for_doy[doy]
          layers_for_doy[doy] += 1
        else
          layers_for_doy[doy] = 1
        end
      end
    end
  end
  (START_DOY..END_DOY).each do |doy|
    # print "doy #{doy}.."; $stdout.flush
    averages[doy] = cumulate_array

    for lati_index in (0..MAX_Y_INDEX)
      for longi_index in (0..MAX_X_INDEX)
        for year in (start_year..end_year)
          # Does this year have this DOY?
          if min = min_grids[year].get_by_index(longi_index,lati_index,doy)
            max = max_grids[year].get_by_index(longi_index,lati_index,doy)
            begin
              averages[doy][lati_index][longi_index] += (max + min) / 2.0
              if lati_index == 2 && longi_index == 22
                # puts "year #{year} doy #{doy}, min #{min} max #{max} avg #{averages[doy][lati_index][longi_index]}"
              end
            rescue Exception => e
              puts "BORKED!\nlati #{lati_index}, longi #{longi_index}, avgs[doy] #{averages[doy].inspect}"
              raise e
            end
          end
        end
      end
    end
  end
  # print "\ndoing the averaging. "
  # Now we have the sum of all the existing daily average temps. Divide by the number
  # of layers that had that DOY.
  (START_DOY..END_DOY).each do |doy|
    "doy #{doy}.."; $stdout.flush
    for lati_index in (0..MAX_Y_INDEX)
      for longi_index in (0..MAX_X_INDEX)
        begin
          raise 'zero in lfd' if layers_for_doy[doy] == 0
          averages[doy][lati_index][longi_index] /= layers_for_doy[doy]
        rescue Exception => e
          puts "BORKED on division! doy #{doy}, lati_index #{lati_index}, longi_index #{longi_index}, lfd[doy] #{layers_for_doy[doy]}"
          puts averages[doy][lati_index].inspect
          raise e
        end
      end
    end
  end
  # puts ""
  averages
end
cell_series(grid_hash,lati_index,longi_index) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 134
def cell_series(grid_hash,lati_index,longi_index)
  series = []
  grid_hash.keys.sort.each do |key|
    series << grid_hash[key][lati_index][longi_index]
  end
  series
end
cumulate(min_temp_hash,max_temp_hash,start_doy,end_doy,base,upper=86) { |min_temp_hash,max_temp_hash,base,upper| ... } click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 77
def cumulate(min_temp_hash,max_temp_hash,start_doy,end_doy,base,upper=86)
  prev = 0
  days_found = 0
  cumulation = (start_doy..end_doy).inject(0) do |sum,doy|
    # Need a better interpolation algorithm here!
    if min_temp_hash[doy] && max_temp_hash[doy]
      days_found += 1
      prev = yield(min_temp_hash[doy],max_temp_hash[doy],base,upper)
    else
      raise "doy #{doy} not found!"
    end
    sum + prev
  end
  [cumulation,days_found]
end
cumulate_array(value=0.0) click to toggle source

return an array so that we can access arr[0..xDim]

# File lib/agwx_biophys/degree_days.rb, line 5
def cumulate_array(value=0.0)
  arr = Array.new(Y_DIM)
  for y in 0..MAX_Y_INDEX
    arr[y] = Array.new(X_DIM,value)
  end
  arr
end
dds_for_year_at_point(p) { |to_fahrenheit(min),to_fahrenheit(max),base,upper| ... } click to toggle source

For a given year, return a hash of DOYs at a grid point with the daily DDs for each, using the passed-in block to calculate a day’s DD

# File lib/agwx_biophys/degree_days.rb, line 254
def dds_for_year_at_point(p)
  calculate_fahrenheit = p[:calculate_fahrenheit] || true
  base = p[:base] || 50.0
  upper = p[:upper] || 86.0
  start_doy = p[:start_doy] || START_DOY # May 1
  [:max_grids,:min_grids,:year,:lati_index,:longi_index].each { |sym| raise "Required param #{sym.to_s} missing" unless p[sym] }
  doy_dds = {} # Will be a Hash keyed by DOY
  
  (start_doy..366).each do |doy|
    min = p[:min_grids][p[:year]].get_by_index(p[:longi_index],p[:lati_index],doy)
    max = p[:max_grids][p[:year]].get_by_index(p[:longi_index],p[:lati_index],doy)
    next unless min && max
    if block_given?
      doy_dds[doy] = yield(to_fahrenheit(min),to_fahrenheit(max),base,upper)
    else
      doy_dds[doy] = modB_DD(to_fahrenheit(min),to_fahrenheit(max),base,upper)
    end
  end
  doy_dds
end
frost?(grid,longi_index,lati_index,doy,frost_value) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 98
def frost?(grid,longi_index,lati_index,doy,frost_value)
  val = grid.get_by_index(longi_index,lati_index,doy)
  return false unless val
  return false unless val != grid.mD.badVal
  return (val <= frost_value)
end
frost_doy(grid,longi_index,lati_index,frost_value=0.0) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 275
def frost_doy(grid,longi_index,lati_index,frost_value=0.0)
  (START_OF_JULY..366).each do |doy|
    min = grid.get_by_index(longi_index,lati_index,doy)
    if min && min <= frost_value
      return doy
    end
  end
  nil
end
frost_doys(min_grids,start_year=START_YEAR,end_year=END_YEAR,frost_value=0.0) click to toggle source

Find the date of first frost – defined as the first day after July 1 with a min temp below 0 C – for each grid point for each year. Return a hash keyed by year, of X by Y arrays of DOYs.

# File lib/agwx_biophys/degree_days.rb, line 107
def frost_doys(min_grids,start_year=START_YEAR,end_year=END_YEAR,frost_value=0.0)
  frost_doys = {}
  for year in (start_year..end_year)
    frost_doys[year] = cumulate_array(NO_DOY)
    # We start with 0 grid cells at frost (it's July, after all). When the number reaches
    # NUM_GRID_CELLS, we know we've seen frost everywhere on the grid and can safely
    # terminate the loop.
    num_frosty_grid_cells = 0
    # Can't use START_DOY or END_DOY here since frost might extend outside the growing season
    (START_OF_JULY..366).each do |doy|
      if min_grids[year].get_by_index(0,0,doy)
        for lati_index in (0..MAX_Y_INDEX)
          for longi_index in (0..MAX_X_INDEX)
            if frost_doys[year][lati_index][longi_index] == NO_DOY && frost?(min_grids[year],longi_index,lati_index,doy,frost_value)
              frost_doys[year][lati_index][longi_index] = doy
              num_frosty_grid_cells += 1
              break if num_frosty_grid_cells >= NUM_GRID_CELLS
            end
            break if num_frosty_grid_cells >= NUM_GRID_CELLS
          end
        end
      end
    end
  end
  frost_doys
end
median(array) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 142
def median(array)
  return nil unless array && array.size > 0 && (array = array.compact).size > 0
  sorted = array.sort
  len = sorted.length
  return (sorted[(len - 1) / 2] + sorted[len / 2]) / 2.0
end
median_pre_frost_remaining_dds(p) click to toggle source

Return a hash, by doy, of the median remaining-DD value for that doy over all years for a point

# File lib/agwx_biophys/degree_days.rb, line 315
def median_pre_frost_remaining_dds(p)
  # hash by year of per-day remaining DDs
  remaining = {}
  [:min_grids,:max_grids,:years,:longi_index,:lati_index].each { |sym| raise "Required param #{sym.to_s} missing" unless p[sym] }
  # Assemble a dataset of remaining DDs for this point; it'll be a hash (by year) of hashes (by doy) of remaining DD
  p[:years].each do |year|
    p[:year] = year # Hack the param hash so it's good to go for pre_frost_remaining_dds_at_point
    remaining[year],frost_doy = pre_frost_remaining_dds_at_point(p)
  end
  # Iterate over each possible DOY, find the median for that DOY. Will pointlessly access many nonexistent DOYs after
  # first frost, but that means we don't terminate prematurely if a DOY happens to be missing
  medians = {}
  (START_DOY..366).each do |doy|
    dds_for_doy = (p[:years].collect {|year| remaining[year][doy]}).compact
    next unless dds_for_doy.size > 0
    medians[doy] = median(dds_for_doy)
  end
  medians
end
modB_DD(min,max,base=50,upper=86) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 28
def modB_DD(min,max,base=50,upper=86)
  min = [base,min].max
  max = [base,max].max
  
  min = [min,upper].min
  max = [max,upper].min
  rect_DD(min,max,base)
end
pre_frost_remaining_dds_at_point(p,proc=nil) click to toggle source

Return a hash of remaining-DD values for a grid point up to the first frost for that year. Most of the params get passed unchanged to dds_for_year_at_point, which handles defaults and such. FIXME: Extend interface to include a passed-in block for DD calcs. For now, just uses dds_for_year_at_point‘s, which is ModB / 50.0 / 86.0

# File lib/agwx_biophys/degree_days.rb, line 300
def pre_frost_remaining_dds_at_point(p,proc=nil)
  
  frost_value = p[:frost_value] || 0.0
  [:max_grids,:min_grids,:year,:lati_index,:longi_index].each { |sym| raise "Required param #{sym.to_s} missing" unless p[sym] }
  fd = frost_doy(p[:min_grids][p[:year]],p[:longi_index],p[:lati_index],frost_value)
  return nil unless fd
  if block_given?
    dds = dds_for_year_at_point(p,&proc)
  else
    dds = dds_for_year_at_point(p)
  end
  [remaining_dds_for_year_at_point(dds,fd),fd]
end
print_samples(yearly_DD_accums,start_year=START_YEAR,end_year=END_YEAR) click to toggle source
rect_DD(min,max,base=50) click to toggle source

Calculate rect DD from min and max for this day. Everything is in Fahrenheit.

# File lib/agwx_biophys/degree_days.rb, line 16
def rect_DD(min,max,base=50)
  rect = ((max + min) / 2.0) - base
  # Never return a negative number
  # puts "rect_dd for #{min},#{max},#{base},#{upper} returning #{rect}"
  rect >= 0.0 ? rect : 0.0
end
rect_DD_from_avg(avg,base=50) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 23
def rect_DD_from_avg(avg,base=50)
  rect = avg - base
  rect >= 0.0 ? rect : 0.0
end
remaining_dds_for_year_at_point(dds_at_point,frost_doy,start_doy=120) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 285
def remaining_dds_for_year_at_point(dds_at_point,frost_doy,start_doy=120)
  remaining_dds = {}
  sum = 0
  frost_doy.to_i.downto(start_doy) do |doy|
    next unless dds_at_point[doy]
    sum += dds_at_point[doy]
    remaining_dds[doy] = sum
  end
  remaining_dds
end
sine_DD(min,max,base=50,upper=DD_MAX_TEMP) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 40
def sine_DD(min,max,base=50,upper=DD_MAX_TEMP)
  alpha = nil
  o1 = o2 = nil
  dd = nil
  avg = (min + max) / 2.0
  #
  # Source Degree-Days: The Calculation and Use
  # of Heat Units in Pest management
  # Original returns an error if min > max
    #           if (TMin > TMax) {
    # RegErr("CalcDayDD:SineDD:","","Tmin > Tmax");
    # return -9999999;
  return upper - base if (min >= upper)
  return 0 if (max <= base)
  return avg - base if (max <= upper && min >= base)

  alpha = (max-min)/2;
  if (max <= upper && min < base)
    o1    = Math.asin( (base-avg)/alpha)
    return M_1_PI*( (avg-base) * (Math::PI / 2-o1) + alpha*Math.cos(o1) )
  end
  
  if (max >  upper && min >= base)
    o2    = Math.asin( (upper-avg)/alpha);
    return M_1_PI*( (avg-base) * (o2+M_1_PI) +
      (upper-base) * (Math::PI / 2-o2) - alpha*Math.cos(o2) )
  end
  
  if (max >  upper && min < base)
    o1    = Math.asin( (base-avg)/alpha);
    o2    = Math.asin( (upper-avg)/alpha);
    return M_1_PI*( (avg-base)*(o2-o1) + alpha*(Math.cos(o1)-Math.cos(o2))
        + (upper-base)*(Math::PI / 2-o2) )
  end
end
to_fahrenheit(celsius) click to toggle source
# File lib/agwx_biophys/degree_days.rb, line 94
def to_fahrenheit(celsius)
  (celsius * (9.0 / 5.0)) + 32.0
end
year_dd_grid(max_grids,min_grids,year) click to toggle source

def cumulate(min_grids,max_grids,averages,start_year=START_YEAR,end_year=END_YEAR,base=10.0,upper=30.0,fahrenheit=false)

# construct a set of DD arrays, one per year
yearly_DD_accums = {}
# print "cumulating. "
for year in start_year..end_year
  yearly_DD_accums[year] = cumulate_array
  # print "#{year}.."; $stdout.flush
  for doy in (START_DOY..END_DOY)
    for lati_index in (0..MAX_Y_INDEX)
      for longi_index in (0..MAX_X_INDEX)
        begin
          min = min_grids[year].get_by_index(longi_index,lati_index,doy) || averages[doy][lati_index][longi_index]
          max = max_grids[year].get_by_index(longi_index,lati_index,doy) || averages[doy][lati_index][longi_index]
          if fahrenheit
            min = to_fahrenheit(min)
            max = to_fahrenheit(max)
            base = to_fahrenheit(base) if base == 10.0
            upper = to_fahrenheit(upper) if upper == 30.0
          end
          if block_given?
            dd = yield(min,max,base,upper)
          else
            dd = rect_DD(min,max,base,upper)
          end
          yearly_DD_accums[year][lati_index][longi_index] += dd
        rescue Exception => e
          puts "cumulate: problem at #{year}, #{doy}, #{lati_index}, #{longi_index}"
          raise e
        end
      end
    end
  end
end
# puts ""
yearly_DD_accums

end

# File lib/agwx_biophys/degree_days.rb, line 249
def year_dd_grid(max_grids,min_grids,year)
  return {1 => 4, 6 => 5}
end