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