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