-
Notifications
You must be signed in to change notification settings - Fork 0
/
gen_potveg_CESM.ncl
187 lines (153 loc) · 5.03 KB
/
gen_potveg_CESM.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
load "~/lib_abrahao.ncl"
; Generate a potential vegetation file in the reference grid using the high resolution map derived using the CESM tools in:
; https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/lnd/clm2/rawdata/dynlnduse_code.050415/
begin
inputfolder = "/home/gabriel/transicao/doutorado/gen_lu_cesm/input/"
outputfolder = "/home/gabriel/transicao/doutorado/gen_lu_cesm/out_potveg/"
reffname = inputfolder + "surfdata.pftdyn_0.9x1.25_simyr1850-2005_c091008.nc"
vegfname = inputfolder + "mksrf_pft_potv_CN.nc"
outfname = outputfolder + "mksrf_pft_potv_CN_0.9x1.25.nc"
; Flip the longitude of the veg map?
flipveg = True
; PFT code range
minpft = 0
maxpft = 16
npft = maxpft - minpft + 1
reffile = addfile(reffname,"r")
vegfile = addfile(vegfname,"r")
; Put proper coordinate variables
vegvar = vegfile->PCT_PFT
vegvar&lat = vegfile->LAT
vegvar&lon = vegfile->LON
; Flip the longitude of the veg file if needed
if (flipveg) then
dumvegvar = vegvar
delete(vegvar)
vegvar = lonFlip(dumvegvar)
delete(dumvegvar)
end if
; Get coordinate variables with coordinate indexing capabilities
veglat = vegvar&lat
veglon = vegvar&lon
veglat!0 = "lat"
veglon!0 = "lon"
veglat&lat = veglat
veglon&lon = veglon
; Reference variable, for the grid
refvar = reffile->PCT_PFT(0,0,:,:)
refvar = refvar@_FillValue
; 2D vars with the coordinates of the boundaries
lats = reffile->LATS
latn = reffile->LATN
lone = reffile->LONE
lonw = reffile->LONW
lonvar = (lone(0,:) + lonw(0,:))/2.0d0
latvar = (lats(:,0) + latn(:,0))/2.0d0
printVarSummary(vegvar)
printVarSummary(refvar)
dims = dimsizes(refvar)
nlatref = dims(0)
nlonref = dims(1)
; Assume all pixels in the subpixel have the same width/height
veglonres = (veglon(dimsizes(veglon)-1) - veglon(0))/(dimsizes(veglon)-1)
veglonw = veglon - (veglonres/2.0)
veglone = veglon + (veglonres/2.0)
veglonw!0 = "lon"
veglone!0 = "lon"
veglonw&lon = veglonw
veglone&lon = veglone
veglatres = (veglat(dimsizes(veglat)-1) - veglat(0))/(dimsizes(veglat)-1)
veglats = veglat - (veglatres/2.0)
veglatn = veglat + (veglatres/2.0)
veglats!0 = "lat"
veglatn!0 = "lat"
veglats&lat = veglats
veglatn&lat = veglatn
; Indexes for reference
indveglatn = ispan(0,dimsizes(veglatn)-1,1)
indveglats = ispan(0,dimsizes(veglats)-1,1)
indveglonw = ispan(0,dimsizes(veglonw)-1,1)
indveglone = ispan(0,dimsizes(veglone)-1,1)
copy_VarCoords(veglatn,indveglatn)
copy_VarCoords(veglats,indveglats)
copy_VarCoords(veglonw,indveglonw)
copy_VarCoords(veglone,indveglone)
; Create the output (pft, lat, lon) variable with PFT fractions
outvar = new((/npft, nlatref, nlonref/), typeof(refvar))
copy_VarCoords_1(vegvar,outvar) ;Copy PFT definition
copy_VarAtts(vegvar,outvar) ;Copy attributes
outvar@comment = "Fraction refers only to nonmissing pixels in the original map, and as such should add up to 100 where there is some land and 0 where theres none"
outvar!1 = "lat"
outvar!2 = "lon"
outvar&lat = latvar
outvar&lon = lonvar
printVarSummary(outvar)
;print(min(lats))
;print(min(latn))
;
;do i = 0, nlatref-1
; print(i + " of " + (nlatref-1))
; do j = 0, nlonref-1
i = 91
j = 239
;Here we get the coordinates of the cells boundaries
latnpix = veglatn({lats(i,j):latn(i,j)})
latspix = veglats({lats(i,j):latn(i,j)})
lonepix = veglone({lonw(i,j):lone(i,j)})
lonwpix = veglonw({lonw(i,j):lone(i,j)})
fulllatnpix = array_append_record(latnpix,veglatn(indveglats({latspix(dimsizes(latspix)-1)})),0)
fulllatspix = array_append_record(veglats(indveglatn({latnpix(0)})),latspix,0)
fulllonepix = array_append_record(lonepix,veglone(indveglonw({lonwpix(dimsizes(lonwpix)-1)})),0)
fulllonwpix = array_append_record(veglonw(indveglone({lonepix(0)})),lonwpix,0)
delete(latspix)
delete(latnpix)
delete(lonwpix)
delete(lonepix)
latspix = fulllatspix
latnpix = fulllatnpix
lonwpix = fulllonwpix
lonepix = fulllonepix
delete(fulllatspix)
delete(fulllatnpix)
delete(fulllonwpix)
delete(fulllonepix)
;print("lats(i,j) = " + lats(i,j))
;print("latn(i,j) = " + latn(i,j))
;print("lonw(i,j) = " + lonw(i,j))
;print("lone(i,j) = " + lone(i,j))
;print(latnpix)
;print(latspix)
;print(lonepix)
;print(lonwpix)
exit
vegpix = vegvar(:,{lats(i,j):latn(i,j)},{lonw(i,j):lone(i,j)})
print(dimsizes(vegpix))
print(vegpix)
print(dim_sum_n_Wrap(vegpix,0))
k = 0
; if (all(ismissing(vegpix))) then
; outvar(:,i,j) = 0.0
; else
; do k = minpft,maxpft
; ;k = 1
; ;print(num(vegpix.eq.k))
; ;print(num(.not.ismissing(vegpix)))
; ;print( 100.0*(todouble(num(vegpix.eq.k)) / todouble( num(.not.ismissing(vegpix))) ))
; outvar({k},i,j) = (100.0*(todouble(num(vegpix.eq.k)) / todouble( num(.not.ismissing(vegpix))) ))
; end do ;k, pfts
; end if ; all(ismissing(vegpix))
exit
delete(vegpix) ; Pixel size can vary
; end do ;j, lons
;end do ;i, lats
; Write output
system("rm " + outfname)
outfile = addfile(outfname,"c")
outfile->vegtype = outvar
;; This is to check that, in the CESM file, PFT fractions dont add up to 100
;poi = reffile->PCT_PFT((/0,155/),:,:,:)
;printVarSummary(poi)
;soma = dim_sum_n_Wrap(poi,1)
;printVarSummary(soma)
;dum_write(soma,"soma")
end