remapcon limitations
Hi all,
I need to find the truth about the limitations of the remapcon command. I have abused this command to generate the area-weighted mean over polygons. In this thread [https://code.mpimet.mpg.de/boards/1/topics/12461?r=12463#message-12463] Uwe stated that for remapcon (or, more properly, the YAC coupler code) to give proper results the grid cells must be strictly convex. The attached image illustrates the nature of my polygons, which are not convex. I need to know 100% whether remapcon results for polygons such as these are incorrect, along with a guess of how incorrect they might be. The results have already been published in an internal technical report and a paper is pending.
The initial project was so successful that I have now been tasked with repeating my calculations ( 
) for many more such polygons and these new polygons are even more complex that the ones in my figure. I realize I may instead have to use the selregion command for my future calculations. The beauty of remapcon is that, ideally, it was weighting and normalizing the model results over the same geographic area. I was using this to trim ocean model grid cell along the shelf break, to get values more appropriate to the continental shelf.
EDIT: text from the YAC paper
"The clipping of two polygons computes the intersection
between them. Typical clipping algorithms assume that all
edges of the polygons are straight lines. Due to the type
of edges supported by YAC, this is not the case here. We
have tried to approximate the latitude circle edges using
sections of great circles, but this increases computation time and
makes the cell concave even though in latitude–longitude
space it is convex. Currently, we use a modified version
of the Sutherland–Hodgman clipping algorithm (Sutherland
and Hodgman, 1974). This algorithm requires that one cell is
either a convex cell that has only great-circle edges or a
rectangular cell consisting of latitude and longitude circle edges.
The second cell can be either convex or concave and can have
any combination of edge types."
My interpretation of this is that my results should be okay. My polygons are not convex, but the ocean model grids cells are.
Thank-you!
Replies (9)
RE: remapcon limitations - Added by Ralf Mueller over 1 year ago
Hi again!
I had a long debugging session with Moritz and our analysis of the current 2.1.0 release version (yac did not change in the last releases AFAIK) is this:
- the clipping algorithm does a triangulation based on the given cell center of your region
- for convex cells this works
- non-convex cells, which have the property that every connection from the cell center each of the corners is part of the cell also work
- completely arbitrary cell shapes do not work, because the triangulation leads to triangles that to not fully overlap with the original cell
Hence, you should get errors with shapes that are critical. possible with something like this:
We also checked the cell area computation, which is used to normalize the individual weights wrt. to the source cells. This computation does take into account, that some triangles might partly lay outside the original cell. So it should be ok for critial shapes.
So much for what the code should be able to do ...
Finally we tested all 4 regions with a constant regular grided field like this:
cdo -info -remapcon,cdo_scrip_grid_ess.nc -const,1.0,global_0.1
In case the weight computation does not work for your regions, some areas should have been counted multiple times, which should lead to value != 1.0. This final values includes weight computation (which should give wrong results) and normalization (which should be ok for your regions). Hence we expected values != 1.0, but we only got 1.0 for all 4 inputs:
% cdo -s -info -remapcon,cdo_scrip_grid_css.nc -const,1.0,global_0.1 [Thu 2022-10-20|14:26:13]
-1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter ID
1 : 0001-01-01 00:00:00 0 1 0 : 1.0000 : -1
% cdo -s -info -remapcon,cdo_scrip_grid_ess.nc -const,1.0,global_0.1 [Thu 2022-10-20|14:26:21]
-1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter ID
1 : 0001-01-01 00:00:00 0 1 0 : 1.0000 : -1
% cdo -s -info -remapcon,cdo_scrip_grid_wss.nc -const,1.0,global_0.1 [Thu 2022-10-20|14:26:33]
-1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter ID
1 : 0001-01-01 00:00:00 0 1 0 : 1.0000 : -1
% cdo -s -info -remapcon,cdo_scrip_grid_gom.nc -const,1.0,global_0.1 [Thu 2022-10-20|14:26:38]
-1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter ID
1 : 0001-01-01 00:00:00 0 1 0 : 1.0000 : -1
To me this a good but surprising results. Of course it strongly depends on the input shapes, which I did not investigate any further.
Possible next steps:
- Create artificial non-convex inputs and check CDOs behaviour on those. This might give a better understanding of the size of the possible error. Instead of remapping the fields we can check the weight computation instead with the
genconoperator. - Check the difference of the results when using the alternative
selregion-based workflow: Since selregion does not involve any interpolation, you might considerable change the data-area compared to your target region. That's why I think it's pretty natural to involve remapping here. But I agree: The errors must be analyzed. - Replace the input region polygons with a (possible unstructured) grid that has only convex cells. Such an input is supposed to work with CDO. You can create such coordinates with libraries like http://gamma.cs.unc.edu/SEIDEL/ or with python's matplotlib like here . Conservative remapping works fine with unstructured grids like this.
Sorry, Brendan, for the inconclusive first answer. I will create an artificial grid like a U and rerun the above tests.
cheers
ralf