clamp: lai/01.read.land_class.split.ncl@0c6405ab2ff4

     1 ;************************************************

     2 ;    Read ascii, Write nc

     3 ;    output: lat: N->S     lon: -180W->180E

     4 ;    class 8 split into class  8 = woody savanna (S. Hem.)

     5 ;                             18 = woody savanna (N. Hem.)

     6 ;    class 8 split into class  9 =       savanna (S. Hem.)

     7 ;                             19 =       savanna (N. Hem.)

     8 ;************************************************

     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

    10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

    11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

    12 ;************************************************

    13 begin

    14 ;---------------------------------------------------

    15 ; final data

    16   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

    17   filo  = "land_class_0.25deg_new.nc"

    18   c = addfile(diro+filo,"c")

    19 ; filedimdef(c,"time",-1,True)

    21   nlat  = 180*4

    22   nlon  = 360*4

    24   t        = new((/nlat,nlon/),integer)

    25   lon      = new((/nlon/),float)

    26   lat      = new((/nlat/),float)

    27 ;---------------------------------------------------

    28 ; input data

    29   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/modis_landcover_qdeg/"

    30   fili = "modis_landcover_class_qd.asc"

    31   b = diri+fili

    33   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")

    34   lat  = (/ lat(::-1) /)                     ; make N->S

    35   lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east")

    36   lon  = (/ lon - 180. /) ; subtract 180 from all values

    37   lon&lon = lon           ; update coordinates

    38   lat&lat = lat           ; update coordinates

    40 ;=============================

    41 ;  name dimensions of t and assign coordinate variables

    42 ;============================

    43   t!0    = "lat"

    44   t!1    = "lon"

    45   t&lat  = lat

    46   t&lon  = lon

    47   t@long_name = "landcover class"

    48   t@units     = ""

    49   t@_FillValue= -9999

    50   t@missing_value= -9999

    52   t = asciiread(b,(/nlat,nlon/),"integer")

    54   nlat = dimsizes(lat)

    55   ny2  = nlat/2

    57   do j = 0,ny2-1

    58      t(j,:) = where((t(j,:) .eq. 8),18,t(j,:))

    59      t(j,:) = where((t(j,:) .eq. 9),19,t(j,:))

    60   end do

    62   c->LAND_CLASS  = t

    63   c->lat  = lat

    64   c->lon  = lon

    66   print (min(t) + "/" + max(t))

    67 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--