source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/examples/test_rmp_esmf/interp_esmf_unstructured_to_unstructured_one_step.ncl @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 4 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 16.2 KB
Line 
1;----------------------------------------------------------------------
2; This script regrids dummy data from a curvilinear grid to an
3; unstructured grid. Both grids are on one NetCDF file.
4;
5; There's an issue with this script if you try to apply SrcMask2D,
6; because it causes the destination grid to be reduced in the number
7; of points (24,572 to 24,067). The fix seems to be in not setting
8; SrcMask2D, but I would like to figure out why this is the case.
9;----------------------------------------------------------------------
10;--------------------------------------------
11;---- Predefined ESMF modules
12;--------------------------------------------
13load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
14load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
15load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
16;
17; Loading a personal ESMF regridding script to provide
18; a needed bug fix for curvilinear to unstructured.
19;
20;load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
21load "./ESMF_regridding.ncl"
22
23begin
24;--------------------------------------------
25;--- PATHS
26;--------------------------------------------
27  RUNDIR="wkdir_esmf_"+str_lower(srcgrid)+"_"+str_lower(tgtgrid)+"_"+str_lower(interp_method)
28  PATH1="/home/evian/coquart/ESMF/ESMF_with_NCL/toy_calcul_weights_with_esmf/RUNDIR/"
29  PATH2="/home/evian/coquart/ESMF/ESMF_with_NCL/toy_calcul_weights_with_esmf/RUNDIR/"
30  PATH1="./"
31  PATH2="./"
32
33; WARNING : DOES NOT WORK YET FOR CONSERVATIVE REMAPPING
34  if (str_lower(interp_method) .eq. str_lower("conserve")) then
35     print("Does not work yet for conservative remapping")
36     exit
37  end if
38;--------------------------------------------
39;--- Source and Target grid Files used with oasis3-mct for SCRIP
40;--------------------------------------------
41  gridsfile = PATH1+"grids.nc"    ; File with grids at OASIS3-MCT description
42  masksfile = PATH1+"masks.nc"    ; File with masks at OASIS3-MCT description
43
44;--------------------------------------------
45;--- Files
46;-------------------------------------------- 
47  funcini     = PATH2+"fldin.nc"     ; File with the initial field
48  funcinterp  = PATH2+"fldou.nc"     ; File with the remapping field, error
49  srcGridName = PATH2+"source_grid_esmf.nc"
50  tgtGridName = PATH2+"target_grid_esmf.nc"
51  wgtFile     = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_" + interp_method + "_withfillValue.nc"
52  wgtFile1    = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_neareststod.nc"
53  wgtFilennei = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_" + interp_method + ".nc"
54  oa3wgtFile  = PATH2+"rmp_esmf_for_oa3mct_"+srcgrid+"to"+tgtgrid+"_" + interp_method + ".nc"
55
56;--------------------------------------------
57;--- Constants
58;--------------------------------------------
59  msk_field_value = -1
60  msk_err_value = -10000
61  weight_fill_val = -50000
62  pi = 3.14159265359
63  dp_length = 1.2*pi
64  dp_conv   = pi/180.
65
66;--------------------------------------------
67;--- Cleaning: remove any pre-existing file
68;--------------------------------------------
69  system("/bin/rm -f "+funcini)
70  system("/bin/rm -f "+funcinterp)
71  system("/bin/rm -f "+srcGridName)
72  system("/bin/rm -f "+tgtGridName)
73  system("/bin/rm -f "+wgtFile) 
74  system("/bin/rm -f "+wgtFile1) 
75  system("/bin/rm -f "+wgtFilennei)
76  system("/bin/rm -f "+oa3wgtFile)
77  system("/bin/rm -f PET0.RegridWeightGen.Log")
78
79;--------------------------------------------
80;--- Coordinate names
81;--------------------------------------------
82  srcgrdlon = srcgrid+"_lon"
83  srcgrdlat = srcgrid+"_lat"
84  srcgrdclo = srcgrid+"_clo"
85  srcgrdcla = srcgrid+"_cla"
86  srcgrdmsk = srcgrid+"_msk"
87  tgtgrdlon = tgtgrid+"_lon"
88  tgtgrdlat = tgtgrid+"_lat"
89  tgtgrdclo = tgtgrid+"_clo"
90  tgtgrdcla = tgtgrid+"_cla"
91  tgtgrdmsk = tgtgrid+"_msk"
92
93;--------------------------------------------
94;--- Read grids, masks, areas arrays
95;--------------------------------------------
96  gfile   = addfile(gridsfile,"r")
97  mfile   = addfile(masksfile,"r")
98 
99  nt=0
100  src_lon  = gfile->$srcgrdlon$(nt,:)
101  src_lat  = gfile->$srcgrdlat$(nt,:)
102  src_clo  = gfile->$srcgrdclo$(:,nt,:)
103  src_cla  = gfile->$srcgrdcla$(:,nt,:)
104  src_msk  = mfile->$srcgrdmsk$(nt,:)
105  tgt_lon  = gfile->$tgtgrdlon$(nt,:)
106  tgt_lat  = gfile->$tgtgrdlat$(nt,:)
107  tgt_clo  = gfile->$tgtgrdclo$(:,nt,:)
108  tgt_cla  = gfile->$tgtgrdcla$(:,nt,:)
109  tgt_msk  = mfile->$tgtgrdmsk$(nt,:)
110
111  nlatsrc=dimsizes(src_lat)
112  nlonsrc=dimsizes(src_lon)
113  nlattgt=dimsizes(tgt_lat)
114  nlontgt=dimsizes(tgt_lon)
115
116;--------------------------------------------
117;--- Calculate the function to regrid
118;--------------------------------------------
119  func_ana_src = 2. - cos(pi*(acos(cos(src_lat*dp_conv)*cos(src_lon*dp_conv))/dp_length))
120
121;--------------------------------------------
122;--- Calculate the function on the destination grid for the error
123;--------------------------------------------
124  func_ana_tgt = 2. - cos(pi*(acos(cos(tgt_lat*dp_conv)*cos(tgt_lon*dp_conv))/dp_length))
125
126;----------------------------------------------------------------------
127;--- Regridding all in one step
128;----------------------------------------------------------------------
129  Opt                   = True
130  Opt@SrcGridName       = srcGridName
131  Opt@SrcGridLat        = src_lat
132  Opt@SrcGridLon        = src_lon
133  Opt@SrcMask2D         = 1-src_msk     ; Don't set this, otherwise you may get error.
134  Opt@DstGridName       = tgtGridName
135  Opt@DstGridLat        = tgt_lat
136  Opt@DstGridLon        = tgt_lon   
137  Opt@DstMask2D         = 1-tgt_msk
138  Opt@DstGridType       = "unstructured"
139
140  Opt@InterpMethod      = str_lower(interp_method)
141  Opt@WgtFileName       = wgtFile
142
143  if(Opt@InterpMethod.eq.str_lower("conserve")) then
144    Opt@SrcGridCornerLat  = ndtooned(src_cla)
145    Opt@SrcGridCornerLon  = ndtooned(src_clo)
146  end if
147
148  if(Opt@InterpMethod.eq.str_lower("conserve")) then
149    Opt@DstGridCornerLat  = ndtooned(tgt_cla)
150    Opt@DstGridCornerLon  = ndtooned(tgt_clo)
151  end if
152
153  Opt@ForceOverwrite    = True
154  Opt@PrintTimings      = True
155  Opt@Debug             = True
156;
157;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
158; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
159  func_regrid = ESMF_regrid(func_ana_src,Opt)
160  printVarSummary(func_regrid)
161  undef_val=func_regrid@_FillValue
162  printVarSummary(undef_val)
163  delete (Opt)
164  func_regrid_msk_fillvalue  = where((1-tgt_msk) .eq. 1, func_regrid, msk_field_value)
165;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
166;--------------------------------------------------------------------------
167;   Do the regridding "neareststod" all in one step to fill blank points
168;--------------------------------------------------------------------------
169  Opt=True
170  Opt@SrcGridName    = srcGridName
171  Opt@DstGridName    = tgtGridName
172  Opt@WgtFileName    = wgtFile1
173  Opt@SrcGridLat     = src_lat
174  Opt@SrcGridLon     = src_lon
175  Opt@DstGridLat     = tgt_lat
176  Opt@DstGridLon     = tgt_lon
177  Opt@DstMask2D      = 1-tgt_msk
178  Opt@SrcMask2D      = 1-src_msk
179
180  Opt@DstGridType       = "unstructured"
181
182  Opt@InterpMethod   = "neareststod"
183;  Opt@Pole           = "none"
184  Opt@ForceOverwrite = True
185  Opt@PrintTimings   = True
186  Opt@Debug          = True
187
188;---You can choose to skip any of these three steps.
189;  Opt@SkipSrcGrid   = True    ; Will assume source grid file has been generated
190;  Opt@SkipDstGrid   = True    ; Will assume target grid file has been generated
191;  Opt@SkipWgtGen    = True    ; Will assume weights file has been generated
192
193  func_neareststod = ESMF_regrid(func_ana_src,Opt)
194  func_neareststod = where((1-tgt_msk) .eq. 1, func_neareststod, msk_field_value)
195;----------------------------------------------------------------------
196;--- Combine interp_method + neareststod
197;----------------------------------------------------------------------
198  wfile          = addfile(wgtFile,"r") 
199;
200  S_regrid       = wfile->$"S"$
201  dst_lat_w      = wfile->$"yc_b"$
202  dst_lon_w      = wfile->$"xc_b"$
203  src_lat_w      = wfile->$"yc_a"$
204  src_lon_w      = wfile->$"xc_a"$
205  dst_address    = wfile->$"row"$
206  src_address    = wfile->$"col"$
207  frac_dst       = wfile->$"frac_b"$
208  numlinks      = dimsizes(S_regrid)
209  dst_size      = dimsizes(dst_lat_w)
210  src_size      = dimsizes(src_lat_w)
211  print("dst_size_w")
212  print(dst_size)
213;
214  wfile1            = addfile(wgtFile1,"r")
215  S_neareststod     = wfile1->$"S"$
216  dst_lat_n         = wfile1->$"yc_b"$
217  dst_lon_n         = wfile1->$"xc_b"$
218  src_lat_n         = wfile1->$"yc_a"$
219  dst_neareststod   = wfile1->$"row"$
220  src_neareststod   = wfile1->$"col"$
221  links_neareststod = dimsizes(S_neareststod)
222  dst_size_n    = dimsizes(dst_lat_n)
223  print("dst_size_n")
224  print(dst_size_n)
225;
226  if (interp_method .eq. str_lower("neareststod")) then
227;--- Just write the new weights into the remapping file
228     print("Nothing to do")
229  else
230      func_oned                = ndtooned(func_regrid_msk_fillvalue)
231      dst_missing              = new((/dst_size/),"integer")
232      dst_missing@_FillValue   = doubletointeger(undef_val)
233
234      count = -1
235      do i=0,dst_size-1
236         if (ismissing(func_oned(i))) then
237            count = count+1
238            dst_missing(count) = i+1
239         end if
240      end do
241; count > 0 : there are some points that need nearest neighbours
242      if (count .ge. 0) then
243      dst_missing_good    = new((/count+1/),"integer")
244      do i=0,count
245         dst_missing_good(i) = dst_missing(i)
246      end do
247      delete(dst_missing)
248      print("count")
249      print(count)
250      print("dst_missing_good")
251;      print(dst_missing_good)
252;
253;--- Replace weights where function = _FillValue by 1. (neareststod)
254      dst_address_new    = new((/numlinks+count+1/),"integer")
255      src_address_new    = new((/numlinks+count+1/),"integer")
256      S_regrid_new       = new((/numlinks+count+1/),double)
257      print("New size numlinks+count")
258      print(numlinks+count+1)
259      dst_address_new(0:numlinks-1)  = dst_address
260      src_address_new(0:numlinks-1)  = src_address
261      S_regrid_new(0:numlinks-1)     = S_regrid
262      do j=0,count
263         do i=0,links_neareststod-1
264            if ( dst_neareststod(i) .eq. dst_missing_good(j) ) then
265                dst_address_new(numlinks+j) = dst_missing_good(j)
266                src_address_new(numlinks+j) = src_neareststod(i)
267                S_regrid_new(numlinks+j)    = S_neareststod(i)
268                break
269            end if
270         end do
271     end do
272     delete(dst_missing_good)
273; count = 0
274     else
275     dst_address_new    = new((/numlinks/),"integer")
276     src_address_new    = new((/numlinks/),"integer")
277     S_regrid_new       = new((/numlinks/),double)
278     dst_address_new    = dst_address
279     src_address_new    = src_address
280     S_regrid_new       = S_regrid
281     end if
282
283     delete(S_regrid)
284     delete(dst_address)
285     delete(src_address)
286
287     igood = ind (.not. ismissing(S_regrid_new))
288     S_regrid = S_regrid_new(igood)
289     src_address = src_address_new(igood)
290     dst_address = dst_address_new(igood)
291
292; endif neareststod
293     delete(S_regrid_new)
294     delete(dst_address_new)
295     delete(src_address_new)
296     end if
297
298;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
299; Recalculate the interpolated function with the new weights with ESMF
300; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
301  if (interp_method .eq. str_lower("neareststod")) then
302;--- Just write the new weights into the remapping file
303     print("Nothing to do")
304  else
305  filennei = addfile(wgtFilennei,"c")
306  globalAtt               = True
307  globalAtt@title         = "ESMF Offline Regridding Weight Generator"
308  globalAtt@normalization = "destarea"
309  if (interp_method .eq. str_lower("conserve")) then
310     globalAtt@map_method = "Conservative remapping"
311     globalAtt@ESMF_regrid_method = "First-order Conservative"
312  end if
313  if (interp_method .eq. str_lower("bilinear")) then
314     globalAtt@map_method = "Bilinear remapping"
315     globalAtt@ESMF_regrid_method = "Bilinear"
316  end if
317  if (interp_method .eq. str_lower("patch")) then
318     globalAtt@map_method = "Bilinear remapping"
319     globalAtt@ESMF_regrid_method = "Higher-order Patch"
320  end if
321  if (interp_method .eq. str_lower("neareststod")) then
322     globalAtt@map_method = "Bilinear remapping"
323     globalAtt@ESMF_regrid_method = "Nearest source to destination"
324  end if
325  globalAtt@conventions   = "NCAR-CSM"
326  globalAtt@domain_a      = "source_grid_file.nc"
327  globalAtt@domain_b      = "destination_grid_file.nc"
328  globalAtt@grid_file_src = "source_grid_file.nc"
329  globalAtt@grid_file_dst = "destination_grid_file.nc"
330  globalAtt@CVS_revision  = "6.3.0r"
331;
332  fileattdef( filennei, globalAtt )
333  vars = getfilevarnames (wfile)
334  nv=dimsizes(vars)
335  print(nv)
336  S_nnei=S_regrid
337  num_links= dimsizes(S_regrid)
338  var_temp=new((/1/),integer)
339  var_temp(0)=1
340  src_address!0="n_s"
341  dst_address!0="n_s"
342  S_nnei!0="n_s"
343  var_temp!0="num_wgts"
344  delete(src_address@_FillValue)
345  delete(dst_address@_FillValue)
346  delete(S_nnei@_FillValue)
347  do i = 0,nv-4
348     filennei->$vars(i)$=wfile->$vars(i)$
349  end do
350  filennei->$"col"$=src_address
351  filennei->$"row"$=dst_address
352  filennei->$"S"$=S_nnei
353  filennei->$"var_temp"$=var_temp
354
355  Opt                = True
356  Opt@Debug          = True
357  Opt@PrintTimings   = True
358  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFilennei,Opt)
359;  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFile,Opt)
360  printVarSummary(func_regrid_new)
361
362  func_regrid_msk_new  = where((1-tgt_msk) .eq. 1, func_regrid_new, msk_field_value)
363
364  delete(var_temp)
365  end if
366;--------------------------------------------
367;--- Calculate the error on the target grid
368;--------------------------------------------
369  func_ana_tgt       = where((1-tgt_msk) .eq. 1, func_ana_tgt, msk_field_value)
370  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
371       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_new)/func_regrid_msk_new))*100, msk_err_value)
372  else
373       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_fillvalue)/func_regrid_msk_fillvalue))*100, msk_err_value)
374  end if
375  print("Max of the error (equal -10000 on land)")
376  printMinMax(error,True)
377  error_min = where(error .eq. msk_err_value, - msk_err_value, error)
378  print("Min of the error (equal 10000 on land)")
379  printMinMax(error_min,True)
380
381;--------------------------------------------
382;--- Output files: to plot ESMF results and create the rmp file for OASIS3-MCT
383;--------------------------------------------
384  ifile = addfile(funcini,"c")
385  ifile->$srcgrdlon$=src_lon
386  ifile->$srcgrdlat$=src_lat
387  ifile->$"field_in"$=(func_ana_src)
388
389  ofile = addfile(funcinterp,"c")
390  ofile->$tgtgrdlon$=tgt_lon
391  ofile->$tgtgrdlat$=tgt_lat
392  ofile->$"field_ana"$=(func_ana_tgt)
393  ofile->$"field_undef"$=(func_regrid)
394  ofile->$"field_msk_undef"$=(func_regrid_msk_fillvalue)
395  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
396     ofile->$"field_ou"$=(func_regrid_msk_new)
397  end if
398  ofile->$"error"$=(error)
399
400;---- Write the rmp file for OASIS3-MCT
401; if bilinear interpolation : num_wgts=1
402; if conservative interpolation : num_wgts=3
403; if nearest neighbour : num_wgts=1
404
405;--- Normalize to obtain Fracarea option for OASIS3-MCT for conservative interpolation
406  if (interp_method .eq. str_lower("conserve")) then
407     do i=0,num_links-1
408        if (frac_dst(dst_address(i)-1) .ne. 0.0) then
409           S_regrid(i)=S_regrid(i)/frac_dst(dst_address(i)-1)
410        end if
411     end do
412  end if
413
414    oa3file = addfile(oa3wgtFile,"c")
415    dst_lat_w!0="dst_grid_size"
416    src_lat_w!0="src_grid_size"
417    src_address!0="num_links"
418    dst_address!0="num_links"
419    oa3file->$"src_grid_center_lat"$=src_lat_w
420    oa3file->$"dst_grid_center_lat"$=dst_lat_w
421    oa3file->$"src_address"$=src_address
422    oa3file->$"dst_address"$=dst_address
423    if (interp_method .eq. "conserve") then
424       var_temp=S_regrid
425       num_links=dimsizes(var_temp)
426       num_wgts=3
427       var = new ((/num_links,num_wgts/),double)
428       var(:,0) = var_temp
429       var(:,1) = 0.
430       var(:,2) = 0.
431       var!0="num_links"
432       var!1="num_wgts"
433       oa3file->$"remap_matrix"$=var
434  else
435       var_temp=S_regrid
436       num_links=dimsizes(var_temp)
437       num_wgts=1
438       var = new ((/num_links,num_wgts/),double)
439       var(:,0) = var_temp
440       var!0="num_links"
441       var!1="num_wgts"
442       oa3file->$"remap_matrix"$=var
443  end if
444
445end
Note: See TracBrowser for help on using the repository browser.