|
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) |
|
20 |
|
21 nlat = 180*4 |
|
22 nlon = 360*4 |
|
23 |
|
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 |
|
32 |
|
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 |
|
39 |
|
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 |
|
51 |
|
52 t = asciiread(b,(/nlat,nlon/),"integer") |
|
53 |
|
54 nlat = dimsizes(lat) |
|
55 ny2 = nlat/2 |
|
56 |
|
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 |
|
61 |
|
62 c->LAND_CLASS = t |
|
63 c->lat = lat |
|
64 c->lon = lon |
|
65 |
|
66 print (min(t) + "/" + max(t)) |
|
67 end |