xESMF
xESMF copied to clipboard
JiaweiZhuang
tripolar grid
I'm trying to regrid the cmip5 data from the ocean model NEMO, which used a tripolar grid (still curvilinear). It seems none of the bilinear, patch and conservative method works. I thought bilinear and patch do not require grid corners but do they require some other specific grid configurations? The conservative method doesn't work because I can't figure out the corners. The corners in the data are specified as 4 vertices but the orientations seem not constant. I mean in the tripolar grid the four faces of the grid are not necessarily oriented east-west or north-south, what is the orientation that is acceptable? The notebook is uploaded here.
https://github.com/JianghuiDu/cmip5
Your problem seems exactly the same as #14.
See https://github.com/JiaweiZhuang/xESMF/issues/14#issuecomment-369010607, https://github.com/JiaweiZhuang/xESMF/issues/14#issuecomment-369686779
@JianghuiDu, I ran into this issue with GFDL data #14 and had to split the domain and regrid each "part", then combine afterwards.
My working example is here: https://github.com/NicWayand/ESIO/blob/master/notebooks/Regrid_GFDL_Forecast.ipynb
And the split function is here: https://github.com/NicWayand/ESIO/blob/master/esio/import_data.py#L278
I have been meaning to write up a more general example in response to #28, hopefully in the next few weeks.
Thanks! This formatting issue is really annoying. It seems most standard ocean outputs are written in lat_b[i,j,v] format where i and j are model grid center indices and v is the index of the vertices. I just tried the NCAR ncl language which by default accepts this format and applies ESMF regridding without any prior treatments.
I just tried the NCAR ncl language which by default accepts this format and applies ESMF regridding without any prior treatments.
Interesting. Do you have NCL examples for this?
They have a few example scripts. Here is one that I modified to regrid the tripolar data using bilinear, conservative and patch. It calls the vertices directly
sos@lat2d = sfile->lat
sos@lon2d = sfile->lon
Opt@SrcGridCornerLat = sfile->lat_bnds ; corners are necessary
Opt@SrcGridCornerLon = sfile->lon_bnds ; for "conserve" method
where the dimension of lat_bnds and lon_bnds is [i,j,4] and the dimension of lat and lon is [i,j]
The data is regridded on to a 1x1 SCRIP grid.
;======================================================================
; ESMF_regrid_6.ncl
;
; Concepts illustrated:
; - Interpolating from one grid to another using ESMF_regrid
; - Interpolating data from a CMIP5 grid to a 1X1 degree rectilinear grid
;======================================================================
; This example is identical to ESMF_all_6.ncl, except it does the
; regridding in one call to "ESMF_regrid". See ESMF_wgts_6.ncl
; for a faster example of regridding using an existing weights file.
;======================================================================
; This example uses the ESMF application "ESMF_RegridWeightGen" to
; generate the weights.
;
; For more information about ESMF:
;
; http://www.earthsystemmodeling.org/
;======================================================================
; This script regrids a CMIP5 grid to a 1.0 degree world grid and
; plots sea water potential temperature on the new grid.
;
; It uses SCRIP for both the CMIP5 and 1.0 degree world grid.
;======================================================================
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;
; This file still has to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
;---Interpolation methods
methods = (/"bilinear","patch","conserve"/)
;---Input file
srcFileName = "sos_Oclim_CNRM-CM5_piControl_r1i1p1_185001-269912-clim.nc"
wgtFile = "CMIP5_2_World_" + methods + ".nc"
;---Get data and lat/lon grid from CMIP5 Grid
sfile = addfile(srcFileName,"r")
sos = sfile->sos(0,:,:)
sos@lat2d = sfile->lat
sos@lon2d = sfile->lon
print(num(ismissing(sos@lat2d)))
print(num(ismissing(sos@lon2d)))
Opt = True
Opt@SrcFileName = "CMIP5_SCRIP.nc" ; source file name
Opt@DstFileName = "World1deg_SCRIP.nc" ; destination file name
Opt@ForceOverwrite = True
Opt@SrcGridCornerLat = sfile->lat_bnds ; corners are necessary
Opt@SrcGridCornerLon = sfile->lon_bnds ; for "conserve" method
Opt@SrcMask2D = where(.not.ismissing(sos),1,0)
Opt@DstGridType = "1x1" ; Destination grid
Opt@DstTitle = "World Grid 1-degree Resolution"
Opt@DstLLCorner = (/-89.75d, 0.00d /)
Opt@DstURCorner = (/ 89.75d, 359.75d /)
;;Opt@PrintTimings = True
;;Opt@Debug = True
;----------------------------------------------------------------------
; Setup for graphics
;----------------------------------------------------------------------
wks = gsn_open_wks("x11","ESMF_regrid") ; send graphics to PNG file
;---Resources to share between both plots
res = True ; Plot modes desired.
res@gsnDraw = False ; Will panel later
res@gsnFrame = False ; Will panel later
res@gsnMaximize = True ; Maximize plot
res@cnFillOn = True ; color plot desired
res@cnFillPalette = "rainbow" ; set color map
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False ; turn off contour labels
res@cnFillMode = "RasterFill" ; turn raster on
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = ispan(30,40,2)
res@mpFillOn = False
; res@trGridType = "TriangularMesh" ; allow missing coordinates
res@gsnAddCyclic = False
res@pmLabelBarWidthF = 0.7
res@lbLabelBarOn = False ; Will do this in panel
res@gsnAddCyclic = False
;---Resources for paneling
pres = True
pres@gsnMaximize = True
pres@gsnPanelLabelBar = True
pres@lbLabelFontHeightF = 0.01
;----------------------------------------------------------------------
; Loop across each method and generate interpolation weights for
; CMIP5 Grid to World Grid
;----------------------------------------------------------------------
plot_regrid = new(dimsizes(methods),graphic)
do i=0,dimsizes(methods)-1
print("Generating interpolation weights from CMIP5 to")
print("World 1 degree grid using the " + methods(i) + " method.")
Opt@WgtFileName = wgtFile(i)
Opt@InterpMethod = methods(i)
;----------------------------------------------------------------------
; Interpolate data from CMIP5 to World 1-degree grid.
;----------------------------------------------------------------------
sos_regrid = ESMF_regrid(sos,Opt)
printVarSummary(sos_regrid)
;----------------------------------------------------------------------
; Plotting section
;----------------------------------------------------------------------
;---Resources for plotting original data
res@tiMainString = "Data on original CMIP5 grid (" + \
str_join(tostring(dimsizes(sos))," x ") + ")"
plot_orig = gsn_csm_contour_map(wks,sos,res)
;---Resources for plotting regridded data
res@tiMainString = "CMIP5 to 1x1-degree grid (" + \
methods(i) + ") (" + \
str_join(tostring(dimsizes(sos_regrid))," x ") + ")"
plot_regrid(i) = gsn_csm_contour_map(wks,sos_regrid,res)
;---Panel two plots
gsn_panel(wks,(/plot_orig,plot_regrid(i)/),(/2,1/),pres)
;---Clean up before next time in loop.
delete(sos_regrid)
end do
end
Interestingly, from NCL docs, NCL only accepts(m, n, 4) but not (m+1, n+1). It essentially follows SCRIP format, not ESMPy's corner representation.
The only place I can find SrcGridCornerLat in NCL source code is the comments in ESMF_regridding.ncl (!), so I am not sure how that conversion is implemented in NCL.
Given that ESMF/ESMPy is able to read SCRIP format, I guess there's some where in ESMF source code that does the conversion (m, n, 4) -> (m+1, n+1). CC @bekozi: do you know where that is & is it possible to use that in Python level?
If ESMF/NCL is essentially doing something as simple as https://github.com/JiaweiZhuang/xESMF/issues/14#issuecomment-369686779, that can be trivially added as an utility function while xESMF's main API can be kept unchanged.
It looks like curvilinear_to_SCRIP and rectilinear_to_SCRIP are the functions that do the reformatting.
This function might be useful. It is what I have traditionally used to do the conversion. It makes some assumptions about the corner positions in the (m, n, 4) array. This is the companion function to revert from ESMF corners. I noticed the docstrings are a little light on these functions so let me know if you have any questions (if the code is relevant).
@bekozi Thanks, that's helpful!
Looks like they are assuming different corner orderings. get_esmf_corners_from_ocgis_corners() assumes:
ref[0, 0] = _corners[ii, jj, 0]
ref[1, 0] = _corners[ii, jj, 1]
ref[1, 1] = _corners[ii, jj, 2]
ref[0, 1] = _corners[ii, jj, 3]
The reverse conversion create_ocgis_corners_from_esmf_corners() assumes:
# Upper left, upper right, lower right, lower left - this is an ideal order. Actual order may differ later.
slices = [(0, 0), (0, 1), (1, 1), (1, 0)]
Which is the default SCRIP convention? From ESMF docs I see:
The grid corner coordinates need to be listed in an order such that the corners are in counterclockwise order.
So [(0, 0), (0, 1), (1, 1), (1, 0)] seems the correct one (counterclockwise), if the indices are expressed in C-ordering (lat, lon). I guess that get_esmf_corners_from_ocgis_corners() assumes Fortran-ordering (lon, lat)?
CF convention, as far as I know, follows from SCRIP which ESMF adheres to closely: http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#cell-boundaries. And, yes, cell boundaries should defined in a CCW manner in ESMF.
I put a short notebook together to demonstrate how the corner conversion functions I linked to are used. This was for my own benefit as well. ocgis corners may be out of compliance with CF...never confirmed it was compliant necessarily (hence the lack of convention in the function names). You may find you need to adjust indexing. In case it's useful, this is where the (m, n, 4) conversion to ESMF corners is used in ocgis-esmpy.
I guess that get_esmf_corners_from_ocgis_corners() assumes Fortran-ordering (lon, lat)?
It does produce Fortran ordering, but it only works off a single coordinate array. Need to call the function twice for lat and lon corners.