Re: Newbie needs help...

[Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>]

"Liam E. Gumley" <Liam.Gumley@ssec.wisc.edu> writes:
> You mentioned that you are using MODIS data at 1 km resolution (1354 x
> 2030 pixels pre granule). One approach we use is to define a global
> equal area grid; we happen to use 25 km x 25 km grid cells. We loop over
> each 1 km pixel, and based on it's lat/lon, we accumulate the following
> statistics in each grid cell:
> 
> Number of observations,
> Sum of observations,
> Minimum observation,
> Maximum observation.
> 
> From these, we can compute mean and standard deviation. it is quite
> straightforward to then resample the equal area grid to an equal angle
> grid which can be visualized in IDL.

Generally speaking I am a scientist who looks *away* from the earth
not towards it, but with satellite-based astronomy there are still
times where one has to worry about where the satellite is.

I am occassionally charged with constructing various instrument
housekeeping data.  Typically we receive multiple measurements at the
same geographic point and I needed to average them.  I wrote this
little routine which constructed the mean and variance of the data.
It's called GEOSPLAT, but there's really nothing "GEO" about it.  It
simply constructs a histogram of LAT/LON data points, weighted by the
measurement values.  Thankfully our satellite is in a low inclination
orbit, so I never really deal with the poles.  This routine just makes
a 2-d histogram regularly sampled in LON/LAT space.

This routine requires a specially modified version of HIST_2D called
HIST_2DR, but don't worry, all I did was add the REVERSE_INDICES
keyword to HIST_2D.  I can forward that if needed.

This application of REVERSE_INDICES is actually not the most
efficient, since it loops through every cell in the output array.
There is another way to do it with far fewer iterations.  If you have
very large image, then one should shift to that.

Craig

function geosplat, lon, lat, val, expo, xrange=xrange0, yrange=yrange0, $
                   xbinsize=xbinsize, ybinsize=ybinsize, minsamp=minsamp, $
                   variance=hh2, exposure=ee

  forward_function arg_present, hist_2dr

  if n_elements(xrange0) LT 2 then xrange0 = [0,359.9999]
  if n_elements(yrange0) LT 2 then yrange0 = [-90,89.9999]
  xrange = xrange0(0:1)
  yrange = yrange0(0:1)
  if n_elements(xbinsize) LT 1 then xbinsize = 1
  if n_elements(ybinsize) LT 1 then ybinsize = 1
  if n_elements(minsamp) LT 1 then minsamp = 1
  if n_elements(expo) LT n_elements(val) then expo = val*0 + 1.
  if arg_present(hh2) then dohh2 = 1 else dohh2 = 0
  
  hh = hist_2dr(lon, lat, reverse=rr0, $
                min1=xrange(0), max1=xrange(1), bin1=xbinsize(0), $
                min2=yrange(0), max2=yrange(1), bin2=ybinsize(0))
  hh = hh * (val(0)*0)
  ee = hh
  if dohh2 then hh2 = hh
  for i = 0, n_elements(hh)-1 do begin
      if rr0(i+1) - rr0(i) GT minsamp(0) then begin
          ex = expo(rr0(rr0(i):rr0(i+1)-1))
          vv = val(rr0(rr0(i):rr0(i+1)-1))
          hh(i) = total(vv*ex)
          ee(i) = total(ex)
          if dohh2 then hh2(i) = total(vv^2*ex)
      endif
  endfor

  hh = hh / (ee>1e-20)
  hh2 = hh2 / (ee>1e-20)
  hh2 = hh2-hh^2

  return, hh
end



-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------