MCD43A3 coordinate info

From: <burakows_at_nyahnyahspammersnyahnyah>
Date: Thu Apr 17 2014 - 11:20:33 MDT

Hi,

I am using the MODIS MCD43A3 BRDF-Adjusted Albedo product and need advice
on retrieving lat/lon coordinate information. I am aggregating 499 daily
files into a single variable to perform time averaging but cannot figure
out how to assign coordinate information to the aggregated variable. I
tried using .he2 as the extension to copyVarCoords, but this did not
appear to work.

I also run out of memory when I convert from short to float using the
scale factor. I am not sure how to overcome this issue to perform
calculations of blue sky albedo from black-sky albedo (BSAsw) and
white-sky albedo (WSAsw). Could I simply wait to apply the scale at the
very end for just my pixels of interest (~230 sites)?

The relevant parts of my script are below. I have attached an MCD43A3
example file. My yellowstone directory with the remaining files is listed
in the script below.

Thanks for your help.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

;============================================================
; The main code
;============================================================

begin

;----------------------------------------------------------------------
; Load MODIS files section
;----------------------------------------------------------------------

;-----MCD43A3 MODIS BRDF-Adjusted Albedo, 2002-2013
        mdir = ("/glade/p/work/burakows/plot/NCL/albedo/MCD43A3/")
        files = systemfunc("cd "+mdir+" ; ls *.hdf")
  print(files)
        nfiles = dimsizes(files)
  print(nfiles)
        f = addfiles(mdir+files,"r")
        dsizes = getfiledimsizes(f[0])
  print(dsizes)
        WSAsw = f[:]->Albedo_WSA_shortwave(:,:)
        BSAsw = f[:]->Albedo_BSA_shortwave(:,:)

;---Open one file as .he2 to get lat/lon info (THIS DID NOT WORK,
fatal:["Execute.c":6321]:variable (Albedo_WSA_shortwave) is not in file
(feos))
        feos = addfile(mdir+files(0)+".he2","r")
        WSAsw_eos = feos->Albedo_WSA_shortwave

        copy_VarCoords(WSAsw_eos,WSAsw)
        copy_VarCoords(WSAsw_eos,BSAsw)

;---Convert from short to float using scale factor (0.001) (THIS DID NOT
WORK, ran out of memory)
        ;WSAsw_f = WSAsw*WSAsw@scale_factor
        ;BSAsw_f = BSAsw*BSAsw@scale_factor

  printVarSummary(WSAsw)

Variable: WSAsw
Type: short
Total Size: 5748480000 bytes
            2874240000 values
Number of Dimensions: 2
Dimensions and sizes: [YDim_MOD_Grid_BRDF | 1197600] x [XDim_MOD_Grid_BRDF
| 2400]
Coordinates:
Number Of Attributes: 10
  long_name : Albedo_WSA_shortwave
  units : albedo, no units
  valid_range : ( 0, 32766 )
  _FillValue : 32767
  scale_factor : 0.001
  add_offset : 0
  scale_factor_err : 0
  add_offset_err : 0
  calibrated_nt : 5
  hdf_name : Albedo_WSA_shortwave

  printVarSummary(BSAsw)

Variable: BSAsw
Type: short
Total Size: 5748480000 bytes
            2874240000 values
Number of Dimensions: 2
Dimensions and sizes: [YDim_MOD_Grid_BRDF | 1197600] x [XDim_MOD_Grid_BRDF
| 2400]
Coordinates:
Number Of Attributes: 10
  long_name : Albedo_BSA_shortwave
  units : albedo, no units
  valid_range : ( 0, 32766 )
  _FillValue : 32767
  scale_factor : 0.001
  add_offset : 0
  scale_factor_err : 0
  add_offset_err : 0
  calibrated_nt : 5
  hdf_name : Albedo_BSA_shortwave

;---Reshape WSAsw and BSAsw to (ndays,2400,2400)
        WSAsw_3D = reshape(WSAsw,(/499,2400,2400/))
        BSAsw_3D = reshape(BSAsw,(/499,2400,2400/))

;----Extract date (YYYYDDD) from MCD43A filename
;MCD43A3.A2013105.h12v04.005.2013126164822.hdf (example filename)
;0123456789012345678901234567890123456789012345 (position)
; 1 2 3 4 (position x10)
     yyyyddd = new(nfiles,integer)
     do n=0,nfiles-1
        tmp_c = stringtochar(files(n))
        yyyyddd(n) = stringtointeger((/tmp_c(9:15)/)) ;YYYYDOY
     end do

  print(yyyyddd)
  printVarSummary(yyyyddd)

;----Convert YYYYDOY to YYYYMMDD
        yyyymmdd = yyyyddd_to_yyyymmdd(yyyyddd)

  print(yyyymmdd)
  printVarSummary(yyyymmdd)

;----------------------------------------------------------------------
; Calculate baseline average monthly albedo by LC type (snow-covered &
snow-free)
; for all "good" quality (QF = 0) section
;----------------------------------------------------------------------

;----Snow-Covered BSA and WSA linear interpolation (performed on unscaled
data)

        ;--- Create new blue-sky file
        BlueSky_SW=new((/nfiles,2400,2400/),float)
        skyl = 0.2

        do i=0,nfiles-1
                BlueSky_SW(i,:,:) = WSAsw_3D(i,:,:)*skyl +
(1-skyl)*BSAsw_3D(i,:,:)
        end do

        BlueSky_SW@_FillValue = 32767
        BlueSky_SW@long_name = "Blue Sky Albedo"
        BlueSky_SW@units = "albedo, unitless"

printVarSummary(BlueSky_SW)

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Apr 17 11:20:44 2014

This archive was generated by hypermail 2.1.8 : Fri Apr 18 2014 - 17:37:58 MDT