歡迎訂閱微信公眾號:『氣象學家』
閱讀建議
讀取和處理了兩種FY-4A和MODIS衛星數據，進行相關產品的繪圖，插值為不同解析度經緯度格點數據並保存為nc格式文件。拋磚引玉，不做更深入的分析，有任何錯誤歡迎大家批評指正。後文附數據腳本獲取方式。
00.前言介紹
01.中國區域FY4_AGRI_L2 OLR原始數據Lambert投影繪圖和腳本
02.中國區域FY4_AGRI_L2 OLR原始數據經緯度網格投影繪圖和腳本
03.中國區域FY4_AGRI_L2 OLR插值數據經緯度網格投影繪圖和腳本
04.東亞區域MOD06_L2 Cloud_Top_Temperature原始數據繪圖和腳本
05.東亞區域MOD06_L2 Cloud_Top_Temperature插值數據繪圖和腳本
06.數據和腳本獲取
07.參考
00.前言介紹工具：NCL、相關地圖底圖包
配料：風雲4A數據FY4_AGRI_L2 OLR產品、MOD06_L2產品
方法：ESMF_regrid等
成品：中國/東亞區域圖、經緯度格點數據（e.g. , 0.1°*0.1°）
FY4A_read_plot_latlon2d_Lambert.ncl
; =============================================================================; Author: Gavin | Affiliation: NJU
; Email : Zhpfu.atm@gmail.com
; Last modified: 2018-06-18 01:15
; Filename: FY4A_read_plot_latlon2d_to_rectilinear_grid.ncl
; Description: Read FY4A full disc data;
; Take OLR for example;
; Mapping to specific area;
; Plotting China area with China
; Choosing different Map Projections;
; To Be Determined(TBD):
; 問題：需要局部微調，經緯度刻度最好四周都有；投影方式暫時只有圓柱等距投影
; Cylindrical Equidistant Projection；
; =============================================================================
; =============================================================================
; Given a start time and a title, this procedure calls get_cpu_time()
; to get the end_time, then prints "elapsed time" information.
; =============================================================================
procedure print_elapsed_time(start_time,title)
local end_time
begin
end_time = get_cpu_time()
print("======================================================================")
print(title + " elapsed time = " + (end_time-start_time) + " CPU seconds.")
print("======================================================================")
end
; =============================================================================
; ===================================MAIN======================================
; =============================================================================
start_time = get_cpu_time()
fpath = "/Users/zhpfu/Downloads/DATA_MAC/00_FY_satellite/"
fname = "FY4A-_AGRI--_N_DISK_1047E_L2-_OLR-_MULT_NOM_20180104000000_20180104001459_4000M_V0001.NC"
f1 = addfile(fpath+fname,"r")
geohdf= "FY4A_OBI_4000M_NOM_LATLON.HDF"
f2 = addfile(fpath+geohdf,"r")
lon2d = f2->Lon(:,:)
lat2d = f2->Lat(:,:)
; Set default value
; 65534 is the area out of Earth
lon2d@_FillValue = 65534
lat2d@_FillValue = 65534
if (any(ismissing(ndtooned(lon2d)))) then
print("Missing longitude coordinates detected")
end if
if (any(ismissing(ndtooned(lat2d)))) then
print("Missing latitude coordinates detected")
end if
; Set the specific area to plot
minlat = 15.
maxlat = 55.
minlon = 72.
maxlon = 136.
olr =short2flt(f1->OLR)
olr@_FillValue = 32766
; Add the attributes
lat2d@units = "degrees_north"
lon2d@units = "degrees_east"
olr@lat2d = lat2d
olr@lon2d = lon2d
olr@units = "W/M2"
olr@coordinates = "lat2d lon2d"
; =============================================================================
; ================================plot=========================================
; =============================================================================
wks_type = "png"
wks_type@wkWidth = 1200
wks_type@wkHeight = 1200
wks = gsn_open_wks(wks_type,"FY4A_OLR_plot_China") ; send graphics to PNG file
gsn_define_colormap(wks, "NCV_jet") ;MPL_jet
; colors = (/ (/1,1,1/),\
; (/0,0,0/),\
; (/1,1,1/),\
; (/0.647,0.953,0.553/),\
; (/0.239,0.725,0.247/), \
; (/0.388,0.722,0.976/), \
; (/0,0,0.996/),\
; (/0.953,0.020,0.933/),\
; (/0.506,0,0.251/)/)
; gsn_define_colormap(wks,colors)
res = True ; plot mods desired
res@gsnMaximize = True ; maximize plot in frame
res@gsnDraw = False
res@gsnFrame = False
res@gsnSpreadColors = True ; Use full color map
res@sfXArray = lon2d
res@sfYArray = lat2d
res@cnFillOn = True ; color fill
res@cnFillMode = "RasterFill" ; Raster mode is much faster
res@cnRasterSmoothingOn = True
; res@cnFillPalette = "gui_default"
res@cnLinesOn = False ; and uses less memory.
res@cnLineLabelsOn = False
res@trGridType = "TriangularMesh" ; Caution!!! can not ignore!
; =============================================================================
; set for map
res@mpLimitMode = "LatLon"
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
res@mpFillOn = True
res@mpDataSetName = "/Users/zhpfu/Downloads/DATA_MAC/FY_satellite/FY4A/NCL-Chinamap/database/Earth..4"
res@mpDataBaseVersion = "MediumRes" ; or "Ncarg4_1"
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"China"/)
res@mpOutlineSpecifiers = (/"China","China:Provinces"/)
res@mpLandFillColor = "white"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@mpFillBoundarySets = "NoBoundaries"
res@mpOutlineBoundarySets = "NoBoundaries"
res@mpNationalLineColor = "black"
res@mpProvincialLineColor = "black"
res@mpGeophysicalLineColor = "black"
res@mpNationalLineThicknessF = 2
res@mpProvincialLineThicknessF = 1
; =============================================================================
; set for the plot
res@cnFillDrawOrder = "PreDraw"
res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF = 100;min(olr) ; set min contour level
res@cnMaxLevelValF = 300;max(olr) ; set max contour level
res@cnLevelSpacingF = 5.0 ; set contour spacing
res@gsnAddCyclic = False
; res@tmXTOn = True
; res@tmYROn = True
res@lbOrientation = "Horizontal" ; vertical labelbar
res@pmTickMarkDisplayMode= "Always"
; res@pmLabelBarOrthogonalPosF = 0.00 ; Move labelbar up
; res@pmLabelBarParallelPosF = 0.00 ; Move labelbar Right
res@pmLabelBarWidthF = 0.65
res@pmLabelBarHeightF = 0.08
res@lbLabelFontHeightF = 0.018
res@lbPerimOn = False
res@lbLabelAutoStride = True ; Clean up labelbar labels.
res@lbBoxLinesOn = False ; No labelbar box lines.
res@lbLabelFontHeightF = 0.01 ; make labels smaller ( default=0.02 )
res@lbBoxEndCapStyle = "TriangleBothEnds" ; set the two-end triangle
res@tiXAxisString = ""
res@tiYAxisString = ""
res@gsnStringFontHeightF = 0.012
res@gsnLeftString = "OLR"
res@gsnRightString = "W/m~S~2~E~";"W/m2"
plot = gsn_csm_contour_map(wks,olr,res) ; create plot
; =============================================================================
; add South China Sea
nhres = res
nhres@gsnMaximize = False
nhres@vpHeightF = 0.18
nhres@vpWidthF = 0.18
nhres@mpMinLatF = 2
nhres@mpMaxLatF = 23
nhres@mpMinLonF = 105
nhres@mpMaxLonF = 123
nhres@lbLabelBarOn = False
nhres@tmXBOn = False
nhres@tmXTOn = False
nhres@tmYLOn = False
nhres@tmYROn = False
nhres@gsnLeftString = ""
nhres@gsnRightString = ""
map_nanhai = gsn_csm_contour_map(wks,olr,nhres)
adres = True
adres@amParallelPosF = 0.499 ; -0.5 is the left edge of the plot.
adres@amOrthogonalPosF = 0.50 ; -0.5 is the top edge of the plot.
adres@amJust = "BottomRight"
plotnh = gsn_add_annotation(plot,map_nanhai,adres)
; add Changjiang and Huanghe river
river = True
river@gsLineThicknessF = 3.0
river@gsLineColor = "blue"
plotrv = gsn_add_shapefile_polylines(wks,plot,"/Users/zhpfu/Downloads/DATA_MAC/00_FY_satellite/FY4A/NCL-Chinamap/cnmap_NetCDF/river.nc",river)
draw(plot)
frame(wks)
print_elapsed_time(start_time,"Plotting Over >>>>>>>")
FY4A_read_plot_latlon2d_to_rectilinear_grid.ncl
04.東亞區域MOD06_L2 Cloud_Top_Temperature原始數據繪圖和腳本MODIS_Full_Area_read_plot_latlon2d_to_rectilinear_grid.ncl
; =============================================================================
; Author: Gavin | Affiliation: NJU
; Email : Zhpfu.atm@gmail.com
; Last modified: 2019-07-19 01:15
; Filename: MODIS_Full_Area_regrid_latlon2d_to_rectilinear_grid.ncl
; Description: Transform xy kilometer to lat lon grid;
; Output data with NetCDF;
; To Be Determined(TBD):
;
; =============================================================================
procedure print_elapsed_time(start_time,title)
local end_time
begin
end_time = get_cpu_time()
print("======================================================================")
print(title + " elapsed time = " + (end_time-start_time) + " CPU seconds.")
print("======================================================================")
end
; =============================================================================
; ===================================MAIN======================================
; =============================================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
start_time = get_cpu_time()
fpath = "/Users/zhpfu/Downloads/DATA_MAC/MODIS/Full_Area/"
fname = "MOD06_L2.A2016095.0340.061.2017326010830.hdf"
f1 = addfile(fpath+fname,"r")
lon2d = f1->Longitude(:,:)
lat2d = f1->Latitude(:,:)
; Set default value
; 65534 is the area out of Earth
lon2d@_FillValue = -999.9
lat2d@_FillValue = -999.9
if (any(ismissing(ndtooned(lon2d)))) then
print("Missing longitude coordinates detected")
end if
if (any(ismissing(ndtooned(lat2d)))) then
print("Missing latitude coordinates detected")
end if
; Set the specific area to plot
minlat = 15
maxlat = 75
minlon = 72
maxlon = 145
; Read cloud temperature.
wv1s = f1->Cloud_Top_Temperature
; Apply scale and offset and convert to double.
CTT0 = wv1s@scale_factor*1.d * (wv1s - wv1s@add_offset)
CTT0@_FillValue = -32768
; Add the attributes
lat2d@units = "degrees_north"
lon2d@units = "degrees_east"
CTT0@lat2d = lat2d
CTT0@lon2d = lon2d
CTT0@units = "K"
CTT0@coordinates = "lat2d lon2d"
print("++++++++++++++++AAAAA++++++++++++++++++")
; Options to set for regridding
interp_method = "bilinear" ; Interpolation method: "bilinear" (default), "patch", or "conserve"
Opt = True
Opt@WgtFileName = "XY_to_rect.nc"
Opt@SrcGridLat = lat2d
Opt@SrcGridLon = lon2d
Opt@SrcRegional = True ; the default
Opt@SrcInputFileName = fpath+fname
Opt@DstRegional = True
Opt@SrcMask2D = where(.not.ismissing(CTT0),1,0) ; Necessary if has
; missing values.
Opt@InterpMethod = interp_method
Opt@DstGridType = "0.05deg" ; destination grid
Opt@DstTitle = "China Grid 0.05 degree resolution"
Opt@DstLLCorner = (/minlat, minlon/) ;;--Change (maybe)
Opt@DstURCorner = (/maxlat, maxlon/) ;;--Change (maybe)
print("++++++++++++++++BBBBB++++++++++++++++++")
Opt@ForceOverwrite = True
Opt@CopyVarCoords = True
Opt@Debug = True
CTT = ESMF_regrid(CTT0,Opt)
printVarSummary(CTT)
print("++++++++++++++++CCCCC++++++++++++++++++")
; ;---Interpolate to a 0.05x0.05 grid
; Opt = True
; Opt@ForceOverwrite = True
; Opt@Debug = True
; Opt@DstGridType = "0.05deg"
; Opt@WgtFileName = wgt_filename
; Opt@SrcMask2D = where(.not.ismissing(CTT),1,0)
; CTT = ESMF_regrid (CTT, opt)
; end if
; l3m_regrid@long_name = CTT@long_name + " (0.05x0.05)"
print("++++++++++++++++DDDDD++++++++++++++++++")
; =============================================================================
; Use the rcm2rgrid_Wrap to regridding
lat = new(1201 , "float",lat2d@_FillValue)
lon = new(1461, "float",lon2d@_FillValue)
lat = ispan(minlat*100,maxlat*100,5)*0.01
lon = ispan(minlon*100,maxlon*100,5)*0.01
lat@_FillValue = -999.9
lon@_FillValue = -999.9
; =============================================================================
lat@units = "degrees_north" ;don't forget to assign the LatLon attributes
lon@units = "degrees_east"
lat!0 = "lat"
lat&lat = lat
lon!0 = "lon"
lon&lon = lon
CTT!0 = "lat"
CTT&lat = lat
CTT!1 = "lon"
CTT&lon = lon
CTT@units = "K"
CTT@coordinates = "latlon"
printVarSummary(CTT)
print("Get subregional data!")
;===================================================================
; Assume variables T, PS and ORO exist and that they have
; associated meta data: (a) coordinate variables time, lev, lat, lon
; and (b) attributes
;===================================================================
nlat = dimsizes(lat)
nlon = dimsizes(lon)
filo = "CTT_China.nc" ; Output file
system("/bin/rm -f " + fpath + filo) ; remove if exists
fout = addfile (fpath + filo, "c") ; open output file
;===================================================================
; explicitly declare file definition mode. Improve efficiency.
;===================================================================
setfileoption(fout,"DefineMode",True)
;===================================================================
; create global attributes of the file
;===================================================================
fAtt = True ; assign file attributes
fAtt@title = "NCL Efficient Approach to netCDF Creation"
fAtt@source_file = ""
fAtt@Conventions = "None"
fAtt@creation_date = systemfunc ("date")
fileattdef( fout, fAtt ) ; copy file attributes
;===================================================================
; predefine the coordinate variables and their dimensionality
; Note: to get an UNLIMITED record dimension, we set the dimensionality
; to -1 (or the actual size) and set the dimension name to True.
;===================================================================
dimNames = (/"lat", "lon"/)
dimSizes = (/nlat, nlon/)
dimUnlim = (/False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim)
;===================================================================
; predefine the the dimensionality of the variables to be written out
;===================================================================
; Here we are using NCL functions to facilitate defining
; each variable's dimension name(s) and type.
; The following could be replaced with explicit, user defined dimension
; names different from those associated with the variable in memory.
; Say, PS(time,lat,lon) in the NCL script. They could be redefined for the file via:
; filevardef(fout, "PS" ,typeof(PS) ,(/"TIME","latitude","longitude"/))
;===================================================================
filevardef(fout, "lat" ,typeof(lat) ,getvardims(lat))
filevardef(fout, "lon" ,typeof(lon) ,getvardims(lon))
filevardef(fout, "CTT" ,typeof(CTT) ,getvardims(CTT))
;===================================================================
; Copy attributes associated with each variable to the file
; All attributes associated with each variable will be copied.
;====================================================================
filevarattdef(fout,"lat" ,lat) ; copy lat attributes
filevarattdef(fout,"lon" ,lon) ; copy lon attributes
filevarattdef(fout,"CTT" ,CTT) ; copy CTT attributes
;===================================================================
; explicitly exit file definition mode. **NOT REQUIRED**
;===================================================================
setfileoption(fout,"DefineMode",False)
;===================================================================
; output only the data values since the dimensionality and such have
; been predefined. The "(/", "/)" syntax tells NCL to only output the
; data values to the predefined locations on the file.
;====================================================================
fout->lat = (/lat/)
fout->lon = (/lon/)
fout->CTT = (/CTT/)
; =============================================================================
; ================================plot=========================================
; =============================================================================
wks_type = "png"
wks_type@wkWidth = 2400
wks_type@wkHeight = 2400
wks = gsn_open_wks(wks_type,"MODIS_CTT_regrid_plot_China") ; send graphics to PNG file
gsn_define_colormap(wks, "NCV_jet") ;MPL_jet
; colors = (/ (/1,1,1/),\
; (/0,0,0/),\
; (/1,1,1/),\
; (/0.647,0.953,0.553/),\
; (/0.239,0.725,0.247/), \
; (/0.388,0.722,0.976/), \
; (/0,0,0.996/),\
; (/0.953,0.020,0.933/),\
; (/0.506,0,0.251/)/)
; gsn_define_colormap(wks,colors)
res = True ; plot mods desired
res@gsnMaximize = True ; maximize plot in frame
res@gsnDraw = False
res@gsnFrame = False
res@gsnSpreadColors = True ; Use full color map
; res@sfXArray = lon2d
; res@sfYArray = lat2d
res@cnFillOn = True ; color fill
res@cnFillMode = "RasterFill" ; Raster mode is much faster
; res@cnRasterSmoothingOn = True
; res@cnFillPalette = "gui_default"
res@cnLinesOn = False ; and uses less memory.
res@cnLineLabelsOn = False
; res@trGridType = "TriangularMesh" ; Caution!!! can not ignore!
; =============================================================================
; set for map
res@mpLimitMode = "LatLon"
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
res@mpFillOn = True
res@mpDataSetName = "/Users/zhpfu/Downloads/DATA_MAC/FY_satellite/FY4A/NCL-Chinamap/database/Earth..4"
res@mpDataBaseVersion = "MediumRes" ; or "Ncarg4_1"
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"China","Japan", "North Korea","South Korea", "Russia"/)
res@mpOutlineSpecifiers = (/"China","China:Provinces"/)
res@mpLandFillColor = "white"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@mpFillBoundarySets = "NoBoundaries"
res@mpOutlineBoundarySets = "NoBoundaries"
res@mpNationalLineColor = "black"
res@mpProvincialLineColor = "black"
res@mpGeophysicalLineColor = "black"
res@mpNationalLineThicknessF = 2
res@mpProvincialLineThicknessF = 1
; =============================================================================
; set for the plot
res@cnFillDrawOrder = "PreDraw"
res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF = 180;min(CTT) ; set min contour level
res@cnMaxLevelValF = 280;max(CTT) ; set max contour level
res@cnLevelSpacingF = 10.0 ; set contour spacing
res@gsnAddCyclic = False
; res@tmXTOn = True
; res@tmYROn = True
res@lbOrientation = "Horizontal" ; vertical labelbar
res@pmTickMarkDisplayMode= "Always"
; res@pmLabelBarOrthogonalPosF = 0.00 ; Move labelbar up
; res@pmLabelBarParallelPosF = 0.00 ; Move labelbar Right
res@pmLabelBarWidthF = 0.65
res@pmLabelBarHeightF = 0.08
res@lbLabelFontHeightF = 0.018
res@lbPerimOn = False
res@lbLabelAutoStride = True ; Clean up labelbar labels.
res@lbBoxLinesOn = False ; No labelbar box lines.
res@lbLabelFontHeightF = 0.01 ; make labels smaller ( default=0.02 )
res@lbBoxEndCapStyle = "TriangleBothEnds" ; set the two-end triangle
res@tiXAxisString = ""
res@tiYAxisString = ""
res@gsnStringFontHeightF = 0.012
res@gsnLeftString = "Cloud_Top_Temperature"
res@gsnRightString = "K";"W/m2"
plot = gsn_csm_contour_map(wks,CTT,res) ; create plot
; =============================================================================
; add South China Sea
nhres = res
nhres@gsnMaximize = False
nhres@vpHeightF = 0.18
nhres@vpWidthF = 0.18
nhres@mpMinLatF = 2
nhres@mpMaxLatF = 23
nhres@mpMinLonF = 105
nhres@mpMaxLonF = 123
nhres@lbLabelBarOn = False
nhres@tmXBOn = False
nhres@tmXTOn = False
nhres@tmYLOn = False
nhres@tmYROn = False
nhres@gsnLeftString = ""
nhres@gsnRightString = ""
map_nanhai = gsn_csm_contour_map(wks,CTT,nhres)
adres = True
adres@amParallelPosF = 0.499 ; -0.5 is the left edge of the plot.
adres@amOrthogonalPosF = 0.50 ; -0.5 is the top edge of the plot.
adres@amJust = "BottomRight"
plotnh = gsn_add_annotation(plot,map_nanhai,adres)
; add Changjiang and Huanghe river
river = True
river@gsLineThicknessF = 3.0
river@gsLineColor = "blue"
plotrv = gsn_add_shapefile_polylines(wks,plot,"/Users/zhpfu/Downloads/DATA_MAC/FY_satellite/FY4A/NCL-Chinamap/cnmap_NetCDF/river.nc",river)
draw(plot)
frame(wks)
print_elapsed_time(start_time,"Game Over >>>>>>>")
插值為0.05°x0.05°格點
06.數據和腳本獲取由於篇幅限制，所有的腳本、測試數據、樣圖，大小約160MB+，公眾號後臺回復關鍵詞獲取：「fy4a」
07.參考1.https://www.ncl.ucar.edu/
2.https://coding.net/u/huangynj/p/NCL-Chinamap/git
3.http://climate2weather.cc
另，中國區域地圖部分，為了表示對貢獻者勞動成果的尊重，若使用該地圖數據繪圖發表論文等，可考慮添加致謝！中文致謝：感謝中國科學院大氣物理研究所黃永傑博士提供的包含正確中國國界和行政區劃的地圖數據
英文致謝：Thank Dr. Yongjie Huang (IAP/CAS) for providing map database (https:
Yong-Jie Huang (IAP/CAS) huangynj@gmail.com
歷史文章推薦（點擊閱讀）
Julia語言在氣候系統模式中的應用
Julia程式語言助力天氣/氣候模式
GRIB格式數據處理
使用Python處理NetCDF格式文件
Nature(2019)-地球系統科學領域的深度學習及其理解
都卜勒雷達反演風場工具-PyDDA| SciPy 2019
並行下載最新ERA-5數據的Python腳本
利用PyCINRAD處理、顯示天氣雷達基數據
Python中如何使用NCL的全部色表Colormaps？
GMT5(The General Mapping Tools)初探の安裝、配置、運行
有任何問題都歡迎交流探討，共同學習進步！
點個試試！ ↓❤↓۞↓➹↓♨↓۞↓