load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" wcStrt = systemfunc("date") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Part1: Reading CSV files(stations_NW.csv, rbscReport_edited5_NW_R24h_Meteomanz.csv) print("Part1: Reading .csv files (stations_NW.csv)") ;(rbscReport_edited5_NW_R24h_Meteomanz.csv) fname_NW = "stations_NW.csv" lines_NW = asciiread(fname_NW,-1,"string") delim = "," ;---Read fields 1-5 from fname_Report station_id_NW = tointeger(str_get_field(lines_NW(1:),1,delim)) icao_code_NW = str_get_field(lines_NW(1:),2,delim) st_name_NW = str_get_field(lines_NW(1:),3,delim) lat_NW = tofloat(str_get_field(lines_NW(1:),4,delim)) lon_NW = tofloat(str_get_field(lines_NW(1:),5,delim)) print(station_id_NW+" "+icao_code_NW+" "+st_name_NW+" "+lat_NW+" "+lon_NW) nSta = dimsizes(station_id_NW) ;---Print the information print("End Part1 Reading .csv files (stations_NW.csv)") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Part2: Interpolating SMAP data for April 2016 to April 2017 print("Part2: Interpolating SMAP files") print("====") diri = "./" fili = systemfunc("cd "+diri+" ; ls SMAP_L3_SM_P_20170*h5") nfili= dimsizes(fili) filesc= str_get_cols(fili,23,30) print(filesc) print("================") pthi = diri + fili print(pthi) print("================") SMAP_fod=new((/nfili,nSta/),float) ;---Open files; 'src_file' is type list src_file = addfiles(pthi, "r") ;---Set group; begin and end with / grp_smrd = "/Soil_Moisture_Retrieval_Data_AM/" var_long_name = "SMAP: Soil Moisture" ; change; default long_name is too long ;---Set variable name varName = "soil_moisture" do nf=0,nfili-1 print("########################################") print(" nf="+nf+": "+fili(nf)) print("########################################") ;---Same lat/lon dimension sizes but maybe different ;---There is no mention _FillValue=-9999.0 in the file. latName = "latitude" lat_path = grp_smrd + latName lat2d := src_file[nf]->$lat_path$ printVarSummary(lat2d) printMinMax(lat2d, 0) print("================") lonName = "longitude" lon_path = grp_smrd + lonName lon2d := src_file[nf]->$lon_path$ printVarSummary(lon2d) printMinMax(lon2d, 0) print("================") ;---There are -9999.0 in the variable ;---Manually add _FillValue lat2d@_FillValue = -9999.0 lon2d@_FillValue = -9999.0 printMinMax(lat2d,0) printMinMax(lon2d,0) print("================") ;---Set variable and group path var_path = grp_smrd + varName var := src_file[nf]->$var_path$ ; curvilinear grid printVarSummary(var) printMinMax(var, 0) print("================") Opt = True ;---Destination file options Opt@DstGridType = "1x1" ; 1x1 degree grid Opt@DstLLCorner = (/-89.75, 0.00 /) Opt@DstURCorner = (/ 89.75, 359.75 /) Opt@GridMask = where(.not.ismissing(var),1,0) Opt@SrcGridLat = lat2d Opt@SrcGridLon = lon2d Opt@InterpMethod = "bilinear" ; default is bilinear ; [ neareststod, nearestdtos ] Opt@ForceOverwrite = True Opt@PrintTimings = True Opt@RemoveSrcFile = True ; remove SCRIP grid source description file Opt@RemoveDstFile = True ; remove SCRIP grid destination " " Opt@RemoveWgtFile = True Opt@NoPETLog = True Opt@Title = "SMAP: "+Opt@InterpMethod+"_to_"+Opt@DstGridType ;+---ESMF regrid wcStrtESMF = systemfunc("date") begTime = get_cpu_time() var_regrid = ESMF_regrid(var,Opt) ; rectilinear grid print("SMAP: ESMF_regrid: CPU:"+Opt@InterpMethod+"_to_"+Opt@DstGridType+" "+(get_cpu_time() - begTime) + " seconds") wallClockElapseTime(wcStrtESMF, "SMAP: ESMF_regrid: Wall:", 0) printVarSummary(var_regrid) printMinMax(var_regrid, 0) print("================") ;---Bilinear interpolation ;---Will happen *ONLY* when target location is surrounded by 4 values/ wcStrt_linint2 = systemfunc("date") begTime = get_cpu_time() SMAP_fo = linint2_points(var_regrid&lon, var_regrid&lat, var_regrid, True, lon_NW, lat_NW, 0) print("SMAP: linint2_points: "+(get_cpu_time() - begTime) + " seconds") wallClockElapseTime(wcStrt_linint2, "SMAP: linint2: Wall:", 0) printVarSummary(SMAP_fo) printMinMax(SMAP_fo, 0) print("================") ;---Extract date ; SMAP_L3_SM_P_20170401_R14010_001.h5 ; 1 2 3 ; 01234567890123456789012345678901234 yyyymmdd = toint( str_get_cols(fili(nf), 13, 20) ) print("************************************************") print(">>> yyyymmdd: "+yyyymmdd+" after ESMF Regrid <<<") print("************************************************") print( sprintf("%6.2f", lat_NW)+" "+sprintf("%6.2f", lon_NW) \ ; array print +" "+sprintf("%6.2f", SMAP_fo)) print("") print("&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&") ;---Plot options begTime = get_cpu_time() pltDir = "./" pltType = "png" pltName = "PLOT."+Opt@InterpMethod+"_to_"+Opt@DstGridType+"."+fili(nf) plot = new(2, "graphic") ;---Plot pltPath = pltDir+pltName wks = gsn_open_wks(pltType,pltPath) res = True ; Plot modes desired. res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnFillOn = True ; color plot desired res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour labels res@cnFillMode = "RasterFill" ; turn raster on res@cnFillPalette = "BlAqGrYeOrReVi200" res@lbLabelBarOn = False ; turn off individual cb's if (varName.eq."soil_moisture") then res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = 0.05 ; set min contour level res@cnMaxLevelValF = 0.95 ; set max contour level res@cnLevelSpacingF = 0.05 ; set contour spacing end if ;---Resources for plotting original (source) data ;---Resources for lat2d, lon2d with _FillValue res@trGridType = "TriangularMesh" res@sfXArray = lon2d res@sfYArray = lat2d ;---Set contour resources res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = 0.05 ; set min contour level res@cnMaxLevelValF = 0.95 ; set max contour level res@cnLevelSpacingF = 0.05 ; set contour spacing ;---Original curvilinear grid res@gsnAddCyclic = False ; True is the default res@gsnLeftString = var_long_name res@gsnCenterString = "Original Grid" plot(0) = gsn_csm_contour_map(wks,var,res) ; original variable ;---Next grids are rectilinear. These resources are not needed. delete([/ res@sfXArray, res@sfYArray, res@trGridType /]) ;---Rectilinear res@gsnAddCyclic = True ; True is the default res@gsnCenterString = "ESMF: "+Opt@InterpMethod+": "+Opt@DstGridType plot(1) = gsn_csm_contour_map(wks,var_regrid, res) ; create plot ; pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@lbLabelFontHeightF = 0.01 ;pres@lbBoxLinesOn = False pres@gsnPanelMainString = fili(nf) gsn_panel(wks,plot,(/2,1/),pres) print("SMAP: Plot: "+(get_cpu_time() - begTime) + " seconds") exit ; for test only one file end do ; end 'nf' print("@@@@@@@@@@@@@@@") wallClockElapseTime(wcStrt, "Processing and Graphics", 0) print("@@@@@@@@@@@@@@@")