In my previous post I’ve described how you can take grid data and turn it into contour data. In this post I’ll show how I implemented this using C#.
The input to my implementation is a simple two-dimensional array of double values (double[][]).
Now I will need a class that will hold the information of a single point on the grid where a contour intersects with a grid line (a ‘hit’ in step II). We’ll call this class IsobarPoint:
public class IsobarPoint {
public int Value { get; set; }
public Coordinate Coordinate { get; set; }
public Point Location { get; set; }
public IsobarDirection Direction { get; set; }
public Isobar Parent { get; set; }
}
As you can see the IsobarPoint contains a field for its value (the threshold values passed), and two types of coordinates –
Coordinate – an [x,y] value of whole numbers, indicating the index of the point on the grid line (for the first hit in the sample image this will be [1,0].
Location – this is an [x,y] value of numbers, indicating the exact location of the point on the area we want to display (i.e. the geo-location). For the first point it might look like [1.75, 0]. This is the data described in Step III.
The IsobarPoint also contains a field for Direction, as described in Step III. It’s values may be East, West, North or South.
The last field points the the Isobar that this IsobarPoint is part of. If it is not set, it means that this point was not processed yet.
Now we need an algorithm which takes the grid data, and creates the IsobarPoints, as described in Steps II and III:
public static IEnumerable CreateIsobars(double[][] data) {
if (data == null) {
return null;
}
var hgrid = new IEnumerable[data.Length][];
var vgrid = new IEnumerable[data.Length - 1][];
for (int x = 0; x < data.Length; x++) {
if (x != data.Length - 1) {
vgrid[x] = new IEnumerable[data[x].Length];
}
hgrid[x] = new IEnumerable[data[x].Length - 1];
for (int y = 0; y < data[x].Length; y++) { var value = data[x][y]; if (x != 0) { vgrid[x - 1][y] = Enumerable.Range(Math.Min((int)value, (int)data[x - 1][y]) + 1, Math.Abs((int)data[x - 1][y] - (int)value)) .Select(v => new IsobarPoint {
Coordinate = new Coordinate(x - 1, y),
Location = new
Point((x - (v - value) / (data[x - 1][y] - value)), y),
Direction = (value > data[x - 1][y]) ? IsobarDirection.East : IsobarDirection.West,
Value = v
}
).ToList();
}
if (y < data[x].Length - 1) { hgrid[x][y] = Enumerable.Range(Math.Min((int)value, (int)data[x][y + 1]) + 1, Math.Abs((int)data[x][y + 1] - (int)value)) .Select(v => new IsobarPoint {
Coordinate = new Coordinate(x, y),
Location = new Point(x, (y + (v - value) / (data[x][y + 1] - value))),
Direction = (value > data[x][y + 1]) ? IsobarDirection.South : IsobarDirection.North,
Value = v
}
).ToList();
}
}
}
return GenerateIsobars(vgrid, hgrid);
}
In this method we’re creating two grids of lists of IsobarPoints – hgrid and vgrid. hgrid will hold all the horizontal IsobarPoints, and vgrid all the vertical ones.
This implementation presumes that the threshold is a whole number, and that the area of the samples start from [0,0] and has a value every 1.0 degree. This can easily be changed (or parameterized).
Notice that hgrid and vgrid are of different sizes from each other, and from the original data grid. This is because grid lines are only between points, so for the hgrid, there is one less column, and for the vgrid there is one less row.
The last line in our method delegates the data constructed in the two grids to another method whose function is to iterate over them, and create the Isobars themselves:
private static IEnumerable GenerateIsobars(IEnumerable[][] vgrid, IEnumerable[][] hgrid) {
// iterate the frame
foreach (var l in vgrid) {
foreach (var i in l[0].Where(i => i.Direction == IsobarDirection.East)) {
yield return new Isobar(i, vgrid, hgrid, false);
}
foreach (var i in l.Last().Where(i => i.Direction == IsobarDirection.West)) {
yield return new Isobar(i, vgrid, hgrid, false);
}
}
foreach (var i in hgrid[0].SelectMany(l => l.Where(i => i.Direction == IsobarDirection.South))) {
yield return new Isobar(i, vgrid, hgrid, false);
}
foreach (var i in hgrid.Last().SelectMany(l => l.Where(i => i.Direction == IsobarDirection.North))) {
yield return new Isobar(i, vgrid, hgrid, false);
}
// find circles
for (int y = 1; y < vgrid.Length - 1; y++) {
for (int x = 1; x < vgrid[y].Length - 1; x++) { foreach (var i in vgrid[y][x].Where(i => i.Parent == null)) {
yield return new Isobar(i, vgrid, hgrid, true);
}
}
}
}
This method iterates over the grids as described in Stage IV (connect the dots), and for each unprocessed IsobarPoint (i.e. has no Parent), it creates a new Isobar. The Isobar class in responsible for tracing the Isobar point, picking up all the points on the resulting path:
public class Isobar {
public Isobar(IsobarPoint first, IEnumerable[][] vgrid, IEnumerable[][] hgrid, bool isClosed) {
Value = first.Value;
IsClosed = isClosed;
Points = GetPoints(first, vgrid, hgrid).ToArray();
}
private IEnumerable GetPoints(IsobarPoint first, IEnumerable[][] vgrid,
IEnumerable[][] hgrid) {
first.Parent = this;
yield return first;
IsobarPoint next, current = first;
while ((next = current.FindNext(vgrid, hgrid)) != null && next != first) {
next.Parent = this;
yield return next;
current = next;
}
}
public IsobarPoint[] Points { get; private set; }
public int Value { get; private set; }
public bool IsClosed { get; private set; }
}
(For those of you which are unfamiliar with C# and the yield return syntax, each yield return denotes an item in the resulting list).
Now all that is missing is the FindNext method in the IsobarPoint class, which can find the next IsobarPoint on the list, as described in Step IV (it is a bit long, but rather repetitive, as it looks for something bit different for each of the possible directions of the IsobarPoint…)
public IsobarPoint FindNext(IEnumerable[][] vgrid, IEnumerable[][] hgrid) {
IEnumerable candidates = null;
switch (Direction) {
case IsobarDirection.East:
if (Coordinate.Y == vgrid[0].Length - 1) {
return null;
}
candidates = (vgrid[Coordinate.X][Coordinate.Y + 1]).Where(x => x.Direction == IsobarDirection.East);
if (Coordinate.Y x.Direction == IsobarDirection.North)));
}
break;
case IsobarDirection.West:
if (Coordinate.Y == 0) {
return null;
}
candidates = (vgrid[Coordinate.X][Coordinate.Y - 1]).Where(x => x.Direction == IsobarDirection.West);
if (Coordinate.X x.Direction == IsobarDirection.South));
}
candidates =
candidates.Concat(
(hgrid[Coordinate.X][Coordinate.Y - 1].Where(x => x.Direction == IsobarDirection.North)));
break;
case IsobarDirection.North:
if (Coordinate.X == 0) {
return null;
}
candidates = hgrid[Coordinate.X - 1][Coordinate.Y].Where(x => x.Direction == IsobarDirection.North);
if (Coordinate.X > 0) {
candidates =
candidates.Concat(
(vgrid[Coordinate.X - 1][Coordinate.Y].Where(
x => x.Direction == IsobarDirection.West)));
if (Coordinate.Y x.Direction == IsobarDirection.East));
}
}
break;
case IsobarDirection.South:
if (Coordinate.X == hgrid.Length - 1) {
return null;
}
candidates = hgrid[Coordinate.X + 1][Coordinate.Y].Where(x => x.Direction == IsobarDirection.South);
if (Coordinate.X x.Direction == IsobarDirection.East));
candidates =
candidates.Concat(
(vgrid[Coordinate.X][Coordinate.Y].Where(
x => x.Direction == IsobarDirection.West)));
}
break;
}
if (candidates == null) {
return null;
}
return candidates
.Where(x => x.Parent == null && x.Value == Value).FirstOrDefault();
}
You can download this code, as well as a WPF window which displays the data according to the stages of the previous post here.
In part III of this post I will discuss some of the math I’ve skipped in part I (for example, the lemmas I used).
Enjoy
Uri Agassi's Blog 
I've been struggling for 2 days to come up with an algorithm to do this, I really appreciate the code & ideas. Thank you
LikeLike
You are very welcome
LikeLike
Great job – save me a LOT of work… if you have a PayPal account I'll happily buy you a beer or six. BTW Any reason for using jagged arrays for your RawData rather than regular 2D arrays? I wonder how easy would it be to include an array of levels as an input?
LikeLike
Thank you very much!
I'm not sure I remember the reason why I chose a jagged array, but it might be over some implementation difficulty I encountered, and decided to go around instead of solve (I don't remember what it was…)
Feel free to take my code and build features on it, and if you produce something cool, give us a link there!
LikeLike
For an alternative approach take a look at: http://paulbourke.net/papers/conrec/
LikeLike
Thanks!
A very impressive paper, seems to me that it might be appropriate for use-cases different than the one I needed to solve, but it produces very nice results 😉
LikeLike
Thank you sir!
I only had to make some minor adjustments to the code (threshold is now parameterized and double instead of int).
I do wonder why the last/first point isn't repeated when the Isobar is closed. When drawing the line, I need to add the missing point manually.
LikeLike
Thank you for your interest!
If you have any improvements to the code you'd like to share I would appreciate it, and update the article in light of it.
LikeLike
Here are my changes to the CreateIsobars function (note that I also changed IsobarPoint so that the Value property is now double). Note that I've only done some basic testing on this but so far so good.
public static IEnumerable CreateIsobars(double[][] data, double stepSize = 1)
{
if (data == null)
{
return null;
}
IEnumerable[][] hgrid = new IEnumerable[data.Length][];
IEnumerable[][] vgrid = new IEnumerable[data.Length – 1][];
for (int x = 0; x < data.Length; x++)
{
if (x != data.Length – 1)
{
vgrid[x] = new IEnumerable[data[x].Length];
}
hgrid[x] = new IEnumerable[data[x].Length – 1];
for (int y = 0; y < data[x].Length; y++)
{
var value = data[x][y];
if (x != 0)
{
double value2 = data[x – 1][y];
int rem0 = (int)(value / stepSize);
int rem1 = (int)(value2 / stepSize);
var newS = new List();
// IF we're crossing a threshold
if (rem0 != rem1)
{
newS = Enumerable.Range(Math.Min(rem0, rem1) + 1, Math.Abs(rem0 – rem1)).Select(a => a * stepSize).ToList();
}
vgrid[x – 1][y] = newS
.Select(v => new IsobarPoint
{
Coordinate = new Coordinate(x – 1, y),
Location = new
Point((x – (v – value) / (data[x – 1][y] – value)), y),
Direction = (value > data[x – 1][y]) ? IsobarDirection.East : IsobarDirection.West,
Value = v
}).ToList();
}
if (y < data[x].Length - 1)
{
double value2 = data[x][y + 1];
int rem0 = (int)(value / stepSize);
int rem1 = (int)(value2 / stepSize);
var newS = new List();
// IF we're crossing a threshold
if (rem0 != rem1)
{
newS = Enumerable.Range(Math.Min(rem0, rem1) + 1, Math.Abs(rem0 – rem1)).Select(a => a * stepSize).ToList();
}
hgrid[x][y] = newS
.Select(v => new IsobarPoint
{
Coordinate = new Coordinate(x, y),
Location = new Point(x, (y + (v – value) / (data[x][y + 1] – value))),
Direction = (value > data[x][y + 1]) ? IsobarDirection.South : IsobarDirection.North,
Value = v
}).ToList();
}
}
}
return GenerateIsobars(vgrid, hgrid);
}
LikeLike
Thanks for your contribution!
LikeLike
@Fat Albert – Thank you, great update 😉
LikeLike
ARGHHH… I hate this comment system… it deleted all that I typed when it logged me in!
Anyway, thanks for this… works great. Mods I made to support GFS Model GRIB2 format where
Lat starts at 90 and decreases to -90, .25 degree increments (721pts)
Lon starts at 0 and increases to 359.75, .25 degree increments (1440pts)
Plot to Google Maps where Lon needs to start at -180 to 180.
In Fat Albert's code the following is update snippets:
public static IEnumerable CreateIsobars(double[][] data, double xResolution = .25, double yResolution = .25, double stepSize = 1)
{
vgrid[x – 1][y] = newS
.Select(v => new IsobarPoint
{
Coordinate = new Coordinate(x – 1, y),
Location = new
Point((x – (v – value) / (data[x – 1][y] – value))* xResolution, y * yResolution),
Direction = (value > data[x – 1][y]) ? IsobarDirection.East : IsobarDirection.West,
Value = v
}).ToList();
hgrid[x][y] = newS
.Select(v => new IsobarPoint
{
Coordinate = new Coordinate(x, y),
Location = new Point(x * xResolution, (y + (v – value) / (data[x][y + 1] – value)) * yResolution),
Direction = (value > data[x][y + 1]) ? IsobarDirection.South : IsobarDirection.North,
Value = v
}).ToList();
In IsobarPoint class, change the Location property:
private Point _location;
public Point Location
{
get { return _location; }
set
{
if (value == null) throw new ArgumentNullException(“value”);
var x = 90 – value.x;
var y = value.y;
if (value.y > 180)
{
y = -(360 – y);
}
_location = new Point(x,y);
}
}
Problem is, I still get horizontal lines for the East/West edge cases as the isoline trys to draw a continuous line to a point on the other side of the map.
I have gone brain dead and could use some help figuring out which parts of the code to change.
Thanks!
LikeLike
Hi Mikla,
Thanks for your interest, and for your contribution!
I'm happy to know my article still helps the community.
From first glance, I believe the problem you are witnessing is because of the new setter in IsobarPoint. Since the algorithm draws from left to right, up down, it actually draws the lines from 0-360 degrees. Since you translate the coordinates to 0-180 and (-180)-0, there is a discontinuity. If you change the coordinates in such a way, you should either _reorder_ the points (-180)-0 to be before 0-180, and then plotting them, or alternatively, plot the map with 0 being the left edge instead of (-180) and plot (-180) to the _right_ of 180 (effectively _not_ translating the Y coordinate, but turning the map instead…)
Uri
LikeLike
After taking a few days, I got things sorted out. I did reorder the points before sending them to the CreateIsobars method. Then in the IsobarPoint I had to shift the data in the Location setter:
private Point _location;
public Point Location
{
get { return _location; }
set
{
if (value == null) throw new ArgumentNullException(“value”);
var x = 90 – value.x;
var y = value.y;
y -= 180;
_location = new Point(x, y);
}
}
Thanks again!
LikeLike
By reading Pressure data from Grib2 File , i have to plot it using (isobars) on google map in my Xamarin.Android App.
Can you plz guide me how to plot isobars on Google Maps V2
LikeLike
Hi Lakhwinder, thanks for reading my post, I hope it helped you. Unfortunately, I'm not familiar with these frameworks, so I'm afraid I would not be much help.
Maybe some other reader might be able to help you…
U.
LikeLike
Hi ,Thanks for the reply .
Your article seems to be very useful and i am sure that it will help me out from my problem .just try to find the way how to use it to plot isobar on google map .
And its not related to any framework . its only related to plotting isobar on Google Map in C#(web desktop , mobile ,any other ).
If you have some idea related to this then that would be a great help for me .
Thanks
LikeLike
Hi , I succeed in making isobar on map by using your Algo, it really saves a lot of time , But now i want the Isobar lines to be smooth not jagged ,so how can i achieve this .
Thanks
LikeLike