Interpolation to grids with <= 3 km grid spacing produces too smooth of height field

[MiCurry]

When interpolating to grids that with 3 km grid spacing, the topo tool produces a height field which is significantly lower resolution than what is expected. I belive this is caused by the double interpolation from the ~1 km GMTED2010 data to the 3 km cube sphere then to the target grid. Then, with ridge analysis the height field is smoothed even further.

Below are two plots of the height field over conus. The first, is the x5.6488066.grid.nc (15-3 km variable resolution mesh) with the 3 km refinement centered over CONUS with the height field that was created with this tool. (I believe the height field outside of the continent lines is due to a lack of landmask).

The second is GMTED2010 dataset interpolated to the x1.65536002.grid.nc (3 km quasi-uniform grid) gird using the MPAS-Model`s init_atmosphere core.

These files can be found on glade at: /glade/scratch/mcurry/cam-mpas-topo-compare

Is there a way interpolate GMTED2010 to a higher resolution cube sphere?

[adamrher]

Hi Miles,

Peter is out this week but I can try to help. I agree the height field is too smooth, and I think it's straightforward to use a higher resolution intermediate cubed-sphere grid if that is indeed the issue. But first, can you confirm some settings in experiment_settings.make? In particular,

ncube_sph_smooth_coarse
nwindow_halfwidth
lregional_refinement

I think I've described to you that you need a refinement factor in the SCRIP file if lregional_refinement=T. And that the other two settings describe the smoothing radius of the base resolution (15km in your case). I'm wondering if it's using the smoothing radius at the base resolution everywhere on the grid?

[MiCurry]

Hey @adamrher thanks for the quick reply. That would be great if you could double check my settings. Here are my experiement settings:

ifeq ($(case),mpas_15-3-conus_Co009_ridge)
  export ncube_sph_smooth_coarse=009
  export output_grid=mpas_15-3
  export grid_descriptor_fname=$(PWD)/cube_to_target/inputdata/grid-descriptor-file/mpasa15-3.conus.desc_SCRIP.210504.nc
  export nwindow_halfwidth=006
  export rdgwin=_Nsw$(nwindow_halfwidth)
  export stitch=-stitch
  export ncube=3000
  export lregional_refinement=.true.
  export rrfac_max=007
  case_found=
endif

You can find my SCRIP file, with rrfac in /glade/work/mcurry/meshes/cam-mpas/mpasa15-3.conus.desc_SCRIP.210504.nc Calculating rrfac for an MPAS grid wasn't exactly as you specified. MPAS does not have a specific dt field for its cells. For instance, it has dcEdge which is the spherical distance between two cells, but not necessarily the cells size. Because MPAS cells are mostly 5, 6, or 7 sided shapes, its hard to calculate accurately determine the grid distance as it might be difference between different vertices between a cell.

However, what I did to calculate rrfac was the calculate the sphere distance between a random grid_corner and the grid_center of the 15-3 km scrip file and then double it, which gives a rough approximation and then divide this calculation by the max distance. An extension to this would be to do this calculation for several corners and take the average.

When looking at the result, the rrfac value was between 1.0 and 7.525. You email said to use 1 for 15 km and 15 / 3 for 3 km.

The code I used to calculate rrfac is here: https://github.com/MiCurry/calc_rrfac/blob/master/calc_rrfac.py

Let me know what your thoughts are on that.

[adamrher]

Those parameters look right. Which is odd b/c your smoothed file (top plot) looks like something you would expect from something closer to 28km grid. As a sanity check, would you be able to plot either the CAM-SE CONUS or the global 14km topo files for comparison?

/glade/p/cesmdata/cseg/inputdata/atm/cam/topo/se/ne30x8_CONUS_nc3000_Co060_Fi001_MulG_PF_RR_Nsw042_c200428.nc
/glade/p/cesmdata/cseg/inputdata/atm/cam/topo/se/ne240pg3_nc3000_Co008_Fi001_PF_nullRR_Nsw005_20171015.nc

As an aside, I also see that for SE, we are setting ncube_sph_smooth_coarse=008 and nwindow_halfwidth=5 for the global 14km grid (ne240pg3), so I'd probably use those parameters for ur mpas-conus grid instead.

I also don't understand how this is a problem with the 1km->3km interpolation. I say that b/c you have the 3km intermediate cubed sphere topofile mapped to the mpas grid, right? And it look more-or-less what we want? If so the 3km source gives the right answer on the mpas grid, and so I suspect it's something in the 3km->target grid smoothing, i.e, cube_to_target.

Lastly, I had thought that the mpas var-res grids were derived from an analytical density function rho(x,y). If so I would think another way to compute the rrfac is to evaluate the density function at cell centers and normalize by the lowest density point on the grid. Maybe this is an outdated workflow ... I recall if from an old mpas paper.

[adamrher]

one more suggestion. As u r trying to just debug this, when running cube_to_target, I would turn off the ridge finder by settting: export rdgwin=_NoAniso ... That should cut down cube_to_target to closer to a half-hour execution time (instead of ~3hrs).

[MiCurry]

Thanks for all the help and suggestions @adamrher.

I just tried changing ncube to be 9000, which would result in a 1 km cube sphere but I'm not sure if this will have any effect as its grabbing the 3km topo-data from: /project/amp/juliob/topo-data/gmted2010_modis-ncube3000-stitch.nc. When calling make, does it only run cube_to_target and not bin_to_cube?

The tool I used to generate the plots above specifically use the MPAS grid field. Do you have scrip files for the files below for plot.nc thats in cube_to_target? I'm not able to locate one for either.

ne30x8_CONUS_nc3000_Co060_Fi001_MulG_PF_RR_Nsw042_c200428.nc
ne240pg3_nc3000_Co008_Fi001_PF_nullRR_Nsw005_20171015.nc

As an aside, I also see that for SE, we are setting ncube_sph_smooth_coarse=008 and nwindow_halfwidth=5 for the global 14km grid (ne240pg3), so I'd probably use those parameters for ur mpas-conus grid instead.

I'll give those a try out.

I also don't understand how this is a problem with the 1km->3km interpolation.

The 1 km -> 3 km interpolation might be correct, but our thinking was that if we are interpolating from a 3 km cube sphere to a 3 km target grid then perhaps we are not 100 % preserving the 3 km details of the cube sphere. i.e. the interpolation of 3 km to any grid spaces would have a resolution of lower than the 3 km cube sphere.

That's at least my interpretation of what might be happening. Perhaps because this tool computes overlap weights that it might be better at preserving detail than I suspect?

I say that b/c you have the 3km intermediate cubed sphere topofile mapped to the mpas grid, right? And it look more-or-less what we want?

My understanding was that the 3km source grid was upon a cube sphere. And cube_to_target takes this cube sphere file and interpolates it to the target grid, regardless of the target_grid. Is that correct?

Lastly, I had thought that the mpas var-res grids were derived from an analytical density function rho(x,y). If so I would think another way to compute the rrfac is to evaluate the density function at cell centers and normalize by the lowest density point on the grid. Maybe this is an outdated workflow ... I recall if from an old mpas paper.

The density function is only used when the variable resolution meshes are created and are not stored within the mesh. I believe they are somewhere on some filesystem, but using them in any means would be a lot of unnecessary work. I don't think we want to be using those.

My suggestion for calculating rrfac is to have this tool calculate it in the way it sees is best. Since the target grid is always a SCRIP file, this tool could use the fields present to calculate this field for every mesh (when lregional_refinement=True) I'd suggest this for two reasons:

1. Reduce user error when calculating or producing it. Calculating the field on the fly will ensure the tool knows exactly what rrfac means
2. By not adding rrfac to the SCRIP file we are leaving the SCRIP file as it is defined and not adding any non-standard fields.

[adamrher]

Hi Miles, I'll respond to you in detail later. But I wanted to first point you to the requested scrip files:
/glade/p/cesmdata/cseg/inputdata/atm/cam/coords/ne0CONUSne30x8_scrip_c200107.nc

the ne240pg3 is not in the repo (@cacraigucar can we do this? You're going to ask me to add metadata aren't ya), but it's here:
/glade/work/aherring/grids/SCRIP_files/ne240pg3_scrip_170628.nc

[adamrher]

I just tried changing ncube to be 9000, which would result in a 1 km cube sphere but I'm not sure if this will have any effect as its grabbing the 3km topo-data from: /project/amp/juliob/topo-data/gmted2010_modis-ncube3000-stitch.nc. When calling make, does it only run cube_to_target and not bin_to_cube?

It does not call bin_to_cube. You need to go into the bin_to_cube folder and generate the executable, and then run the run.sh shell script (this doesn't take long at all). You could experiment with ncube=9000 for debugging purposes, but it's not really a solution b/c 3km was chosen to provide SGH30 for the Beljaars pbl form drag scheme. That is, Beljaars is expecting a subgrid variance for length scales below 3km.

But also note that the "stitch" modifier is important in /project/amp/juliob/topo-data/gmted2010_modis-ncube3000-stitch.nc.. The 1km GMTED dataset has an error over the Antarctic peninsula, and so Julio stitched in a corrected topo for that region into the 3km dataset. We don't have the scripts in the TopoCESM2 to reproduce the corrected 3km dataset. I am in the process of fixing this by stitching in the correction to the 1km dataset ... and eventually I'll include the stitching scripts so that the entire package is reproducible. This also means the current SGH30 over the Antarctic peninsula is wrong, which I'm not sure anyone thought about, but even more motivation to get the fix into the 1km dataset.

The 1 km -> 3 km interpolation might be correct, but our thinking was that if we are interpolating from a 3 km cube sphere to a 3 km target grid then perhaps we are not 100 % preserving the 3 km details of the cube sphere. i.e. the interpolation of 3 km to any grid spaces would have a resolution of lower than the 3 km cube sphere.

That's at least my interpretation of what might be happening. Perhaps because this tool computes overlap weights that it might be better at preserving detail than I suspect?

It might be instructive to layout the workflow in cube_to_target. (1) remap the rrfac to the 3km grid, (2) smooth the topo fields on the 3km grid using the smoothing radius determined by ncube_sph_smooth_coarse, which refers to a smoothing radius in terms of the number 3km grid points. (3) map the smoothed field to the target grid.

So you should expect the fields to be smoother than the 3km grid (although I suppose we would need to ensure that ncube_sph_smooth_coarse >= rrfac_max) [edit for clarity - the smoothing radius is scaled by 1/rrfac_max so that there is less smoothing over the refined region ... basically we are trying to ensure that the smoothing radius is about 1.8dx everywhere, where dx is the grid spacing of the target grid, and obviously dx is smaller in the refined region]. What I think might be worth checking is whether rrfac is taking on values of on the order of rrfac_max over CONUS. If this is not defined properly than the smoothing radius will be too larger in the CONUS region, which what first thought of when I looked at your plots. If you want, you can parse through a version of cube_to_target.F90 that I have (/glade/work/aherring/grids/TopoCESM2/cube_to_target/cube_to_target.F90) where I added in a subroutine to write of the rrfac after its mapped to the 3km grid:

!+++ARH
!
! write netCDF file
!
subroutine wrt_cube(ncube,terr_cube,rrfac_cube,output_file)
  use shr_kind_mod, only: r8 => shr_kind_r8
  use shared_vars
  implicit none

You can search initials ARH to find relevant mods, and how I'm calling it in the script.

[adamrher]

If you want a head start for plotting the 3km grid, you could use my ncl script:
/glade/work/aherring/grids/TopoCESM2/cube_to_target/output/plot_rrfac.ncl

Although its kind of messy and you'll probably have to uncomment some stuff.

[adamrher]

I forgot to respond to your suggestion to compute rrfac as a part of the TopoCESM2 package. Both of your points are good ones. One problem I can think of is that the spectral element grid is a special case. The SCRIP grid only provides the computational grid, which is highly anisotropic. This is because each element contains unequally spaced GLL nodes, and we invent control volumes around those nodal points that results in large differences in control volumes in a single element. Here's an example of what the control volumes look like in a particular coarse uniform spectral element grid:

Because of this, I compute the rrfac on the element grid, which is not available from the SCRIP grid file. The tool that does this is in our toolkit (https://github.com/ESMCI/Community_Mesh_Generation_Toolkit), which is our feeble attempt to standardize this for cam-se grids. But I agree with you that it's not ideal for us to be messing with the well established definition of a scrip grid file.

[MiCurry]

@PeterHjortLauritzen @adamrher thanks again for the meeting today.

Here is the plot of rrfac for the 15-3 grid with its refinement over CONUS:

It uses the code in calc_rrfac and calculates rrfac using the sphere distance from the cell center to a randomly selected cell corner (the first grid corner in the cells list) and then takes the ratio of the max computed distance to each distance.

A few things to note about the plot above that might be noteworthy. First is that the SCRIP grid that we use for topography is on a unit sphere. All fields are scaled to be on a unit sphere. In standalone MPAS, we blow the fields up by an earth radius when we apply static, terrestrial fields (terrain, snow-albedo, vegetation fraction etc.).

The other thing is that all latitude and longitude values are in radians and not degrees. I had to convert the latitude and longitude into degrees. Is it possible that that might have an effect?

I think the strange banding at the top of the refined region (and other places) might be matte pattern caused by the plot, but I'm not sure. However, the several darker spots outside the refine region are consistent with the MPAS grid areas. I'm generating a plot of the MPAS grid areas field for a comparison.

Is it possible that the radian latitude and longitude values have an effect or the fact that the SCRIP grid is on the unit sphere?

[adamrher]

@MiCurry can you point me to the SCRIP grid file and the plot_rrfac.ncl script you used to plot this? I want to try a few things to see if the ugly parts are a plotting artifact.

Is it possible that the radian latitude and longitude values have an effect or the fact that the SCRIP grid is on the unit sphere?

The area in the SCRIP files should be in radians^2, and assumed a sphere of unit radius 1 (so the global integral should be 4 * pi * 1^2). The lat lon values should be in degrees. I'm looking at scrip file that we know works w/ the topo software and I can see that longitudes use the convention [0,360] instead of [-180,180]:
/glade/p/cesmdata/cseg/inputdata/atm/cam/coords/ne0CONUSne30x8_scrip_c200107.nc

Which convention do your converted values correspond to?

[PeterHjortLauritzen]

Thanks Miles and Adam! rrfac looks reasonable. It seems likely that the lat-lon's are 180 degrees off => you might be using rough topography 180 degrees displaced (so not over the US but in the low resolution area and smooth topography in the refinement area).

[MiCurry]

@adamrher I'm using the scrip file found in:/glade/scratch/mcurry/rrfac. The SCRIP I used to make the plot is there as well.

The lat and long values in the scrip file are from [0, 2pi], and translated they are from [0, 360]. Perhaps we can update this tool to translate radians to degrees? We could check the units of each field and translate accordingly. SCRIP files can either have radians or degrees as units for its coordinates:

"The unit of the coordinates can be either "radians" or "degrees"."

The areas are in radians^2.

[PeterHjortLauritzen]

@adamrher: can we overwrite height with refinement factor on the 3km intermediate cubed-sphere so we can map height (refinement factor) to the target grid and see if it is displaced or not?

[adamrher]

@MiCurry the ugliness was a plotting issue. I've found that CellFill doesn't like cylindrical equidistant plots for some reason. If I change it to an orthographic projection it's much happier:

      res@mpProjection = "Orthographic"
      res@mpCenterLatF =    30.
      res@mpCenterLonF =  -100.