A Fantastic Approach for Easy Contour Drawing and Simple High-Low Labeling

Simplifying techniques for contour locating and high-low centre labeling

Charles Luk (cluk19693@gmail.com)

Update: 25th January 2026 on a simplification by removal of a few redundant function codes in appendix B is pending.

Posting on 7th January 2026 :- Complete revision of the entire presentation in SECTION B on High-Low Labeling

Posting date : 11 November 2025

Owner of the website : https://subtropicalrainstorms.ca/2025/08/13/a-study-on-200hpa-disturbance-excitation-behind-a-12-hr-synoptic-scenario-for-imminent-extreme-rainstorm-detection-in-the-subtropical-region/

The presentation in the paper is largely programming-language-independent, except for the two Appendices, which, as extensions of the subject material, implement a walkthrough of Visual C++ program code compiled with Microsoft Visual Studio 2022.

The document, which offers an alternative approach to contour drawing by providing an algorithm for contour detection, linking, and ridge/trough center labeling, starts by converting a 2.5 ° latitude-longitude input grid to an x-y coordinate grid on a Mercator or Polar Stereographic projection. The next step subdivides the coarse row-column grid into a fine mesh grid through carrying out the following steps: 1) divide grid rows into sub-rows called raster, and down each column from the first to the last, interpolate from column-at-row intercepts into column-at-raster intercepts, and 2) divide grid columns into sub-columns, and across each raster from the top to the bottom, interpolate from raster-at-column intercepts into raster-at-sub-column intercepts. Dividing each fine mesh grid value by a contour interval at each grid point, a contour-interval-number integer field is obtained.

Introduction

In the nineteen-eighties, the contour linking method, which plays the crucial role in contour drawing, was based on a raster-cutting algorithm. By examining the grid point values on a raster, the contour-interval-number changes at successive columns are obtained.  Each change location column is at a point of contour intersection on the raster. By arranging intersection locations from left to right on each raster before linking them pairwise from an upper raster to the next lower raster, contour segments between each raster pair are obtained. This method will not be discussed further here.

SECTION A. Contour Drawing

The method discussed was inspired by a University of Toronto Computer Science assignment on simulating water flow from an underground tunnel entrance to its exit. The simulation was generalized here to water flow which gushes from a hole on the ground anywhere inside an empty tunnel (with any number of entrances or exits), any irregularly shaped cave, or closed cavity (with no entrance or exit), each bounded by the wall, eventually reaches all areas of the continuous space.

In the first place, a swath of rasters is generated from the fine-mesh contour-interval-number integer field as a 2D floor plan of the cave. If at any grid (I,J)  the contour-integer-number differs from that of its neighbor grid, it is on a contour which simulates a wall; otherwise, the grid (I,J) represents a white space on the floor plan. The 4 sides of the swathe rectangle represent the extent to which searching (simulating water flow) stops.

Simulation is performed by looping through a fine-mesh contour-interval-number integer field in the swath. At each(the primary) coordinate (I,J) where the grid value is not 255(which means the grid is unvisited), a global variable is set equal to the primary coordinate (I,J)’s value as the comparison target value for the search through the white space. A function call is made, passing the primary coordinate (I,J) as the parent parameter. Inside the function, the grid value at parameter (I,J) is compared against the global variable. If unequal (changing interval number), coordinate (I,J) is a contour intersection and should be appended to a contour intersection queue with immediate return to its call parent. If equal (remarks: inside the first level call by primary parent, the answer is always yes since it is a self comparison), it suggests coordinate (I,J) is in the continuous white space as its parent, and a function call to each neighbor coordinates (except the parent) is made by passing the its (I,J) as the test-location parameter and its own (I,J) as the referring-parent parameter. Because the function call becomes recursive, the result can generate an indefinite number of branches until every individual function returns on 1) encountering a contour intersection, 2) hitting the swathe boundary, or 3) arriving at a visited coordinate (I,J) whose value equals 255. In any case, all return calls made to the parent eventually converge to the primary parent.

Following is a screen(1 and 1/2) shot copy of the recursive function:

Function call made at each primary coordinate (I,J) will result in the detection and storage of contour coordinates in a queue. By linking these coordinates in proper sequence, all the contour lines drawn, which may or may not reach the top, left, right, or bottom of the swathe boundary, would enclose a continuous white space inside which the primary coordinate (I,J) lies. These contours may consist of contours with single or double contour-interval-integer values. Since contour intersections are encountered in a haphazard order during the search, the coordinates of different contour segments of equal or different contour-interval-integer values may occur in the queue in random order. If only reordering the segments without sorting intersection coordinates in the queue, the plotting result may end up in error, as shown in the next screenshot. First, broken or dotted contour lines would appear throughout the chart. Secondly, linking 2 contour intersections in the reverse direction, from right to left, would cause the contour label to be printed upside down, with an upward offset from the gap (refer to the top-left-hand corner of the chart).

Below with sorting applied to the same case (which is tricky), the chart looks good.

SECTION: B. Maxima(High) and Minima(Low) Labeling

All centres of maxima (high or ridge inclusive) and minima (low or trough inclusive) can be located on the chart by simple comparison with surrounding field values. But in most cases, this would lead to duplicating max or min labels. The following shows two examples of duplicating labels.

First example MSL Pressure

Second example 300 hPa contour height

Apart from contour drawing, the recursive function approach for searching contour intersections in Section A can be applied to High/Low labeling. An integer grid field on contour interval number over a basic contour interval number is copied from the master grid onto the grid m_map. Thus, the relationship of the m_map field’s grid value to the actual contour value equals

Actual contour value = (field’s grid value + basic contour interval number) X contour interval width

The array members of Highs are handled top-down. That means High array’s subscribe I loops from total#high to 2. After converting the H(i)’s contour(say MSL pressure) to its contour interval #, it is set equal to a global integer m_examGrade. The Recursive function is called by passing the H(i)’s row and column coordinates(i,j) as the test location parameters. The last argument ‘extype’ is set to +1 for H( as in this case) or -1 for L. Inside the function, the grid content at the test location parameter is compared with m_examGrade. If equal, the location content is set equal to 255. and the same function calls are made at the neighboring locations (coordinates). The succeeding locations will repeat the same process until a chart boundary is reached or a location’s content is not equal to m_examGrade. The function then relinquishes further neighbor calls and eventually backs up to the first parent. The result is a continuous grid domain with a grid value of 255. For each H(i), an inner High array with subscribe k loops from i-1 to 2. if a High array’s H(k)’s coordinates (i,j) on m_map is tested equal to 255, a duplication of H(i) occurring over the same domain is found, resulting the cancellation of either H(i) or H(k) from the list. The result of patching H(i) on m_map is then copied to grid m_patchmap for storing all individual High’s( and later Low’s) patch results from m_map for printing. For C++ implementation, please refer Appendix B.

Function CopyPatch():

First example of corrected single H/L Labeling MSL Pressure

Second example of corrected single H/L Labeling 300 hPa contour height

But this is not the end of story on label printing. Error in label printing does not only confine to duplicating labels over a single domain. The are 2 other common cases of mis-labeling. First, a col is an area or zone surrounded by high or low pressure areas. Therefore it shouldn’t be designated as a high pressure or low pressure area and hence should be label-free. Secondly very often, a low centre label which has no individual enclosing contour at the trailing end of anther deeper low is enclosed by the deeper low’s outer contour forming a complex low system. The same case occurs with some high centre label. This is superfluous as shown in the following example.

It has to be established that each label must be enclosed either by one closed contour or by broken single-value contour fragments with 1(up to 3) open endings on the chart edge boundary. Therefore, a label inside a col or complex H/L does not meet the single-value contour requirement. Any patch formed by calling the ApplyPatch() function on such a label should not be accepted in the subsequent CopyPatch() function call. But there is a simple solution. By introducing a global boolean variable ‘m_bPatchWarn’ in the GradeMap class and use it in the ApplyPatch() function, ‘m_bPatchWarn’ can be set equal to true for any single-value-contour violation in the course of building the patch. Inside the CopyPatch() function call, ‘m_bPatchWarn’ can be used to validate the resulting patch.

The left example shows Highs labelled in the cols over the East China Sea and over Korea, an a complex Low over southwest China are removed. The right example shows the complex Lows across the western Pacific, the philippines, and the South China Sea are reduced to a simple Low.

Summing up, the section above shows that displaying a single label for each high/low/trough/ridge centre, and removing the label from the trailing end of complex H’s/L’s and the col from the chart, could improve chart interpretation.

Appendix A: VC++ code for SECTION A on Contour Drawing

Summing up, upon return from the primary function Apply_b_Patch call at each grid(I,J), the following work will be completed for the white space area enclosing the primary grid: contour intersection queue sorting, queue elements’ contour-interval-number re-grouping, contour intersection retrieving/linking, and turning point labeling. Upon completing the loop over all grids (I,J), drawing all contour lines with turning point labels within the swathe area would be complete. Detailed implementation using Visual C++ code is shown in the screenshots that follow.

VC++ code for Appendix A for contour drawing ends here.

Listing of the full program continues with the start of High-Low Labeling:

Appendix B: VC++ code in High, Low, Trough & Ridge Labeling

As mentioned in SECTION B, all Maximum and Minimum use the synonyms High and Low, respectively, in the discussion. All processing on the High array below applies to the Low array.

Before printing each H/L label, the content of the H/L’s location coordinates(col,row) on the patchmap is tested for -255 as assigned in the function CopyPatch() which is displayed in Section B. All unwanted H/L labels location with content equal to 255 will be removed from the print queue by the following codes.

if (*(Gmp->m_patchmap + Gmp->m_col * (row-1) + col)== 255{/*test made on the print queue*/                      
m_LowRow.RemoveAt(i);
m_LowCol.RemoveAt(i);
}

Listing for high-low labeling

VC++ code for Low Labeling is similar to that of High Lebeling hence the Appendix B for VC++ code listing ends here. The program listing continues to show output-to-printer code.:

Acknowlegement

I would like to dedicate my study to Met O 12(c) of the U.K. Meteorological Office, who provided me with an attachment training on this subject in the early nineteen-eighties under Mr. Adrian Todkill and Mr. P. Cockrell. They were the forerunners of contour chart makers. Without the groundbreaking contribution of their fundamental development in their contour model, my endeavor could not have been realized.

An intriguing Approach for Contour Drawing and Highs, Lows, Troughs, and Ridges Single Labeling

Implementation of new techniques for contour locating and for eliminating duplicate high/low labels in the subtropical regions

Charles Luk (cluk19693@gmail.com)

Another major update (revision) to Appendix B, with a new project title, is pending for the next posting, scheduled for later this month.

The major update(re-written) to Appendix B was completed on 9 December 2025.

Last update 30/11/2025:- 3 ‘copymap’ functions added to the end of Appendix B

Posting date : 11 November 2025

Owner of the website : https://subtropicalrainstorms.ca/2025/08/13/a-study-on-200hpa-disturbance-excitation-behind-a-12-hr-synoptic-scenario-for-imminent-extreme-rainstorm-detection-in-the-subtropical-region/

The presentation in the paper is largely programming-language-independent. However, the two Appendices, which are extensions of the subject material, consist of the techniques implementation through a brief walkthrough of Visual C++ program code compiled on the Microsoft Visual Studio 2022 platform.

The document, which offers an alternative approach to contour drawing by providing an algorithm for contour detection, linking, and ridge/trough center labeling, starts with converting a two-and-a-half-degree latitude-longitude input grid to an x-y coordinate grid on a Mercator or Polar Stereographic projection. The next step subdivides the coarse row-column grid into a fine mesh grid through carrying out the following steps : 1) divide grid rows into sub-rows called raster, and down each column from the first to the last, interpolate from column-at-row intercepts into column-at-raster intercepts , and 2) divide grid columns into sub-columns, and across each raster from the top to the bottom, interpolate from raster-at-column intercepts into raster-at-sub-column intercepts . Dividing each fine mesh grid value by a contour interval at each grid point, a contour-interval-number integer field is obtained.

Introduction

In the nineteen-eighties, the contour linking method, which plays the crucial role in contour drawing, was based on a raster-cutting algorithm. By examining the grid point values on a raster, the contour-interval-number changes at successive columns are obtained.  Each change location column is at a point of contour intersection on the raster. By arranging intersection locations from left to right on each raster before linking them pairwise from an upper raster to the next lower raster, contour segments between each raster pair are obtained. This method will not be discussed further here.

SECTION A. Contour Drawing

The method discussed was inspired by a University of Toronto Computer Science assignment on simulating water flow from an underground tunnel entrance to its exit. The simulation was generalized here to water flow which gushes from a hole on the ground anywhere inside an empty tunnel (with any number of entrances or exits), any irregularly shaped cave, or closed cavity (with no entrance or exit), each bounded by the wall, eventually reaches all areas of the continuous space .

In the first place, a swath of rasters is generated from the fine-mesh contour-interval-number integer field as a 2D floor plan of the cave. If at any grid (I,J)  the contour-integer-number differs from that of its neighbor grid, it is on a contour which simulates a wall; otherwise, the grid (I,J) represents a white space on the floor plan. The 4 sides of the swathe rectangle represent the extent to which searching (simulating water flow) stops.

Simulation is performed by looping through a fine-mesh contour-interval-number integer field in the swath. At each(the primary) coordinate (I,J) where the grid value is not 255(which means the grid is unvisited), a global variable is set equal to the primary coordinate (I,J)’s value as the comparison target value for the search through the white space. A function call is made, passing the primary coordinate (I,J) as the parent parameter. Inside the function, the grid value at parameter (I,J) is compared against the global variable. If unequal (changing interval number), coordinate (I,J) is a contour intersection and should be appended to a contour intersection queue with immediate return to its call parent. If equal (remarks: inside the first level call by primary parent, the answer is always yes since it is a self comparison), it suggests coordinate (I,J) is in the continuous white space as its parent, and a function call to each neighbor coordinates (except the parent) is made by passing the its (I,J) as the test-location parameter and its own (I,J) as the referring-parent parameter. Because the function call becomes recursive, the result can generate an indefinite number of branches until every individual function returns on 1) encountering a contour intersection, 2) hitting the swathe boundary, or 3) arriving at a visited coordinate (I,J) whose value equals 255. In any case, all return calls made to the parent eventually converge to the primary parent.

Following is a screen(1 and 1/2) shot copy of the recursive function:

Function call made at each primary coordinate (I,J) will result in the detection and storage of contour coordinates in a queue. By linking these coordinates in proper sequence, all the contour lines drawn, which may or may not reach the top, left, right, or bottom of the swathe boundary, would enclose a continuous white space inside which the primary coordinate (I,J) lies. These contours may consist of contours with single or double contour-interval-integer values. Since contour intersections are encountered in a haphazard order during the search, the coordinates of different contour segments of equal or different contour-interval-integer values may occur in the queue in random order. If only reordering the segments without sorting intersection coordinates in the queue, the plotting result may end up in error, as shown in the next screenshot. First, broken or dotted contour lines would appear throughout the chart. Secondly, linking 2 contour intersections in the reverse direction, from right to left, would cause the contour label to be printed upside down with an upward off-the-gap displacement (refer to the top left-hand corner of the chart).

Below with sorting applied to the same case (which is tricky), the chart looks good.

All centres of maxima(high or ridge inclusive) and minima(low or trough inclusive) can be located on the chart by simple comparison with surrounding field values. But in most cases, this would lead to duplicating max or min labels. The following shows two examples of duplicating labels.

First example MSL Pressure

Second example 300 hPa contour height

SECTION: B. Maxima or Minima Single Labeling

Apart from contour drawing, the recursive function approach for searching contour intersections in Section A can be modified to eliminate duplicate max(in high or ridge), min(in low or trough) centre labels into a single label. A screenshot of such a function is shown as follows.

All discusions below use the term High for Maximum and Low for Minimum. The array members of Highs are handled top-down. That means High array’s subscribe I loops from total#high to 2. After setting the global variable m_examGrade equal to the value at the High(I)’s coordinates (I,J) on the fine mesh contour-interval-integer grid field (hereafter called a patchmap), the recursive function ApplyPatch call is made from the High(I)’s coordinates (I,J) as the primary parent. Upon returning from the primary call, all cells (I,J) with value m_examGrade on the patchmap before the call would be set to 255, defining a domain area enclosing High(I). Then the successive High at (I-1) to (1) on this patchmap is tested for equality to 255. If true, it means the high centre is a duplicate in the same domain. Before elimination by setting its location centre to 0, the high should swap with the top high if it is closer to the domain centre. After decrementing the top high’s subscript by 1, the process will be repeated until the top high’s subscript reaches the bottom position. The foregoing discussion on High array processing is applied to Low array processing except the results of processing are entered by two different function calls with Gmp->Copy254high() for labeling Maxima and Gmp->Copy252low() for labeling Minima. For how this is implemented in C++, please refer to Appendix B.

Following shows the result of single center labelling:

First example MSL Pressure

Second example 300 hPa contour height

This shows that displaying a single label for high/low/trough/ridge centre on the chart could improve chart interpretation. Hence, it seems desirable to do single labeling.

Appendix A: Contour Drawing

Summing up, upon return from the primary function Apply_b_Patch call at each grid(I,J), the following work will be completed for the white space area enclosing the primary grid: contour intersection queue sorting, queue elements’ contour-interval-number re-grouping, contour intersection retrieving/linking, and turning point labeling. Upon completing the loop over all grids (I,J), drawing all contour lines with turning point labels within the swathe area would be complete. Detailed implementation using Visual C++ program codes is shown in screenshots follows.

Appendix B: High Low Trough & Ridge Labelling

As mentioned in SECTION B, all Maximum and Minimum use the synonyms High and Low, respectively, in the discussion. All processing on the High array applies to the Low array; however, the output is handled separately by two main function calls: Gmp->Copy254high() is used specifically to label Maxima, while Gmp->Copy252low() is used to label Minima. These functions assign the processed results to the appropriate context for further analysis.

First call CopyMap(true) to create and init a low-resolution m_patchmap for storing all individual High’s patch results from m_map. For each High, a loop over the ‘High’ pixel array begins with the previously calculated coordinate on the high-resolution (fine-mesh) grid, which is then converted to the low-resolution grid on m_map. Loop through the array High(i) indexed from total#high to 2. It begins with copying m_master to m_map on each H(i), followed by the primary function call ApplyPatch(High(i)’s col, High(i)’s row, 0, 0) from parent High(i). The purpose of this call is to set all cells connected with High(i) inside the same domain of contour interval(grade) as High(i), to 255 on m_map. It also determines the patch rectangle’s 4 boundaries enclosing the 255 cells. The resulting m_map is then subject to a function call FillZero(bound ‘on the side or bottom edge’, extype ‘set to 1 for high or -ve 1 for low’) to solve the problem later on double labeling a pair of adjacent ridge/trough on both sides of their common boundary in the subtropical zone/region. Thus it is necessary to prepare m_map by calling a recursive function PopulateZero() inside FillZero() which sets, all remaining unpatched cells (i.e. not equal to 255) within the patch rectangle on the m_map, to 254 (for cells with High’s adjacent grade), 252(for cells with Low’s adjacent grade), or 0(all other contour grades). How the solution of double labeling can be realised is described in the function Gmp->Copy252low() below, near the end of the Appendix.

A second loop iterates over the High(k) array, with k indexed from i-1 to 2, to check whether the coordinates (I,J) of H(k) on the patched map are set to 255. This check uses the patched map to determine if H(k) already shares a domain with H(i). If so, it means H(k) duplicates H(i) within the same area, so one of them—whichever is farther from the domain center—is removed, and the preferred High is moved to the current i position. This deduplication step ensures accurate mapping of significant maxima.

When the k-indexed loop is complete, call Gmp->Copy254high() or Gmp->Copy252low(). This step moves H(i)’s or L(i)’s m_map, accumulating results on m_patchmap. In Copy252low(), entries in m_map conflicting with double labeling under adjacent trough/ridge pairs in subtropical regions have their sign reversed from -ve to +ve. This ensures they are removed before printing, as follows:

           if (*(Gmp->m_patchmap + Gmp->m_col * row + col)>0){ /*test made on the print queue*/                      
m_LowRow.RemoveAt(i);
            m_LowCol.RemoveAt(i);
            }

Fill Zero()

Populate zero()

Looping through High Array

}}}

Acknowlegement

I would like to thank Met O 12(c) of the U.K. Meteorological Office, which provided me with an attachment training on this subject in the early nineteen-eighties under Mr. Adrian Todkill and Mr. P. Cockrell. Upon a foundation modeled on an early version of CMAP developed by MetO 12, I built my first contour program.