lai/01.read.land_class.ascii.ncl
changeset 0 0c6405ab2ff4
equal deleted inserted replaced
-1:000000000000 0:44b415ab3437
       
     1 ;************************************************
       
     2 ;    Read ascii, Write nc
       
     3 ;    output: lat: N->S     lon: -180W->180E                         
       
     4 ;************************************************
       
     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
       
     8 ;************************************************
       
     9 begin
       
    10 ;---------------------------------------------------
       
    11 ; final data
       
    12   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
       
    13   filo  = "land_class_0.25deg.nc"
       
    14   c = addfile(diro+filo,"c")
       
    15 ; filedimdef(c,"time",-1,True) 
       
    16        
       
    17   nlat  = 180*4 
       
    18   nlon  = 360*4
       
    19   
       
    20   t        = new((/nlat,nlon/),integer)
       
    21   lon      = new((/nlon/),float)
       
    22   lat      = new((/nlat/),float)
       
    23 ;---------------------------------------------------
       
    24 ; input data
       
    25   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/modis_landcover_qdeg/"
       
    26   fili = "modis_landcover_class_qd.asc"                    
       
    27   b = diri+fili
       
    28   
       
    29   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
       
    30   lat  = (/ lat(::-1) /)                     ; make N->S
       
    31   lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east")
       
    32   lon  = (/ lon - 180. /) ; subtract 180 from all values 
       
    33   lon&lon = lon           ; update coordinates
       
    34   lat&lat = lat           ; update coordinates
       
    35 
       
    36 ;=============================
       
    37 ;  name dimensions of t and assign coordinate variables
       
    38 ;============================  
       
    39   t!0    = "lat"
       
    40   t!1    = "lon"
       
    41   t&lat  = lat
       
    42   t&lon  = lon
       
    43   t@long_name = "landcover class"
       
    44   t@units     = ""
       
    45   t@_FillValue= -9999         
       
    46   t@missing_value= -9999         
       
    47 
       
    48   t = asciiread(b,(/nlat,nlon/),"integer")
       
    49 
       
    50   c->LAND_CLASS  = t
       
    51   c->lat  = lat
       
    52   c->lon  = lon
       
    53 
       
    54   print (min(t) + "/" + max(t))
       
    55 end