Lets create some hexagons in #GeometryMeetsData EP5

Right. Where were we. Oh and apologies for the delay: when you have kids at daycare, you have a higher chance of contracting Flu, especially when the days get shorter and everything gets colder.

In the last Episode EP4 I had put everything in a carousel and I did not find it all that practical, as there is not a lot of space to put explanations nor code. Also, I have gotten some feedback that apparently mixing a horizontal title page with verticle followup pages is definately not a good idea on LinkedIN. Thanks for pointing this out, I am new to this, I had no idea :)

In this #GeometryMeetsData episode 5, we will start creating lists of hexagon cordinates. Next up, in Episode 5, we will be finally converting these lists to Well-Known-Text geometry representations (#wkt) and show it in the Icon Map Visual.

• EP4: discuss several ways to generate hexagon coordinates in #PowerQuery
• EP5: We discuss ways to generate a WKT representation in #PowerQuery
• EP5: We show how to plot it on an icon map visual in #PowerBI

Creating hexagon coordinates

Before we begin, we have to assume that we already created a grid with center coordinates for the hexagonal grid that we are going to plot. This will be explained later, but I think it is more 'depict'-ive if we show something on our screen already in the next episode, rather then in 3 episodes from now.

In EP1 on the basics, we showed that we will by default use a pointy top hexagonal shape. It has a certain location (x,y) and size measured from it's center to the top, or also known as the edge_length. In the H3 system, the average edge lengths are prescribed per level, but that is not really relevant here. The hexagon can also be rotated with a certain angle.

What we are looking for are the coordinates of the hexagon corners, as this allows us later to draw vertices between them. Obtaining these coordinates is (in its basis) a matter of simple trigonometrics: multiply the size of the hexagon (the green line) with the cosine of the angle in radians to find x1, and use the sine to find the y1 coordinate. This gives you the location of the first corner in relation to its center, so you need to add the center x and y coordinates respectively.

In PowerQuery code, the above would look like the following

// Convert the angle to radians for trigonometric calculations.
angle_rad = angle_deg * Number.PI / 180,

// Calculate the x and y coordinates of the corner.
corner_x = x + size * Number.Cos(angle_rad),
corner_y = y + size * Number.Sin(angle_rad),         

Now we only need to find the following 5 corners and luckily that is easy: we can add 60 degrees to angle_deg and calculate again. And then add 120, and then all the way until you have touched upon all corners.

In Powerquery there are many ways to achieve this, and I happened to start with one of the more complex ones in my process. I removed all comments that I have in my code to make it more readable in this article, but in the published version (later) all code is commented exhaustively

Hexagon Trial1: List.Generate

Take note that I calculate 7 (!) sets of coordinates for 7 corners in this case. This was done for later: when creating the polygon WKT representation you need to come back to the original corner and close the loop. (You could also add the first corner to the end instead of looping through 7 corners).

My thought was that looping through all 7 corners with a list.generate, starting from corner 0 to corner 6 would work pretty well.. It would allow the 60 degree offset per corner to be implemented well enough, and it would return me a list of records (x,y) where each row in the list represented a corner of the hexagon. But.... It had an iterator. And although it worked, I am not really fond of that.

let
    pointy_hex_corners_V1 = (x as number, y as number, size as number, rotation_degrees as number) =>
        let
            corners = List.Generate(
                () => [i=0],
                each [i] <= 6,
                each [i=[i]+1],
                each 
                    let
                        angle_deg = 60 * [i] - 30 + rotation_degrees,
                        angle_rad = angle_deg * Number.PI / 180,
                        corner_x = x + size * Number.Cos(angle_rad),
                        corner_y = y + size * Number.Sin(angle_rad),
                        result = [x=corner_x, y=corner_y]
                    in
                        result
            )
        in
            corners
in
    pointy_hex_corners_V1
        

Hexagon Trial2: List.Transform

Similarly, you can 'iterate' through a list using List.Transform. This time no complex list.generate, but a simple {0..6} written out for clarity. It does exactly the same, but apparently after timing it a couple of times, already 20%-ish faster than Trial 1.

let
    pointy_hex_corners_V1_Revised = (x as number, y as number, size as number, rotation_degrees as number) =>
        let
            angles = {0, 1, 2, 3, 4, 5, 6},
            corners = List.Transform(
                angles,
                each 
                    let
                        angle_deg = 60 * _ - 30 + rotation_degrees,
                        angle_rad = angle_deg * Number.PI / 180,
                        corner_x = x + size * Number.Cos(angle_rad),
                        corner_y = y + size * Number.Sin(angle_rad)
                    in
                        [x = corner_x, y = corner_y]
            )
        in
            corners
in
    pointy_hex_corners_V1_Revised
        

Hexagon Trial3: Precalculated angles

Now it got more interesting. I had the feeling that the trigonometric calculations, the Number.PI, Number.Cos etc were 'quite' expensive calculations. Not that it was slow (on my pc 10 seconds for 1M hexagons) but it could be a lot faster, or so I thought. And I was right. When replacing the whole trigonometric calculation with a set of precalculated cosine/sine representations the formula sped up like 3-5 times.

let
    pointy_hex_corners_V3 = (x as number, y as number, size as number) =>
        let
            precalculated = {
                [cos = 0.8660254037844387   , sin = -0.5 ],
                [cos = 0.8660254037844387   , sin = 0.5  ],
                [cos = 0                    , sin = 1.0  ],
                [cos = -0.8660254037844387  , sin = 0.5  ],
                [cos = -0.8660254037844386  , sin = -0.5 ],
                [cos = 0                    , sin = -1.0 ],
                [cos = 0.8660254037844387   , sin = -0.5 ]
            },
            corners = List.Transform(
                precalculated,
                each [x = x + size * [cos], y = y + size * [sin]]
            )
        in
            corners
in
    pointy_hex_corners_V3
        

Now, this comes at a cost. We had to give up on our angle parameter. But this is not a problem for our #dataanalysis cause: H3 has a certain fixed angle at a certain location, and between levels it has a delta of about ~19.1 degrees. In EP4 we discussed the initialization, and the Init record initialized there dynamically based on the desired H3 level precalculates these trigonometrics already for us. So we only need to pass the record to the function.

let
    pointy_hex_corners_V4 = (x as number, y as number, init as record) =>
        let
            size = init[hex_size],
            precalc_trigonometrics = init[precalc_pointyhex_cornertrigonometrics],
            corners = List.Transform(
                precalc_trigonometrics,
                each [x = x + size * [cos], y = y + size * [sin]]
            )
        in
            corners
in
    pointy_hex_corners_V4
        

??Melissa de Korte please support me with your expertise here: is it wise to reference/pass an entire record in this case? Or could you better pass the hex_size and precalc_pointyhex_cornertrigonometrics as seperate arguments? My thought was: for now it's easier to do it like this, because it's clear. But what if Init is really big? Is the PQ parser smart enough to realize you only use 2 parts of it here and adjust? Or would this internally work like C/C++ pointers and no data is passed at all so no performance loss? Questions questions...

Hexagon adjustment for the curves of the earth

In EP2 it was revealed that the earth is not flat. I was not threatened by flat-earthers so that went pretty well. Or maybe flat-earthers don't do PowerBI. Who knows!

But that being said, the above generates hexagons perfectly fine on a FLAT surface and not on the earth using latitude longitude, where they get vertically stretched depending on the desired latitude you put them on. Poor things.

To compensate for that, our trusty EP4 discussed initialization comes back in the game. Do note, that in this version i started sending (x,y) coordinates to the function as records instead of seperate variables. The main difference with V4 above, is the correction to degree from meters in latitude and longitude. This fixes everything. In fact, the images in EP2 with the stretched hexagon was created by temporarily turning off the correction within the hexagon and the grid creation.

Note that the function gets the hexagon center in latitude and longitude, already corrected in V5.

let
    pointy_hex_corners_V5 = (hexcenter as record, init as record) =>
        let
            size = init[hex_size],
            hexcenter_lat = hexcenter[latitude],
            hexcenter_lon = hexcenter[longitude],
            precalc_trigonometrics = init[precalc_pointyhex_cornertrigonometrics],
            lat_meter2degree = init[lat_meter2degree],
            lon_meter2degree = init[lon_meter2degree],
            
            corners = List.Transform(
                precalc_trigonometrics,
                each [hexcorner_lat = hexcenter_lat + size * [cos] * lat_meter2degree, hexcorner_lon = hexcenter_lon + size * [sin] * lon_meter2degree]
            )
        in
            corners
in
    pointy_hex_corners_V5
        

The code for the precalc and correction (plus the inverse) is given as follows:

precalc_pointyhex_cornertrigonometrics  = List.Transform(
    {0..5},
    each 
        let
            angle_deg = 60 * _ + grid_angle,
            angle_rad = angle_deg * Number.PI / 180,
            cos_val = Number.Cos(angle_rad),
            sin_val = Number.Sin(angle_rad)
        in
            [cos=cos_val, sin=sin_val]
)

lat_meter2degree = 1 / 1000 / 111.32

lon_meter2degree = 1 / 1000 / (111.32 * Number.Cos(Number.From(center_lat) * Number.PI / 180))

lat_degree2meter = 1 / lat_correction

lon_degree2meter = 1 / lon_correction        

The function in use

The function in its use is really easy once it is defined. A simple call, and there you have the lists of coordinates! In the next episode, we will use these lists to generate the WKT objects.

let
    getHexCorners= (hexcenter as record, init as record) =>
        let
            size = init[hex_size],
            hexcenter_lat = hexcenter[latitude],
            hexcenter_lon = hexcenter[longitude],
            precalc_trigonometrics = init[precalc_pointyhex_cornertrigonometrics],
            lat_meter2degree = init[lat_meter2degree],
            lon_meter2degree = init[lon_meter2degree],
            
            corners = List.Transform(
                precalc_trigonometrics,
                each [hexcorner_lat = hexcenter_lat + size * [cos] * lat_meter2degree, hexcorner_lon = hexcenter_lon + size * [sin] * lon_meter2degree]
            )
        in
            corners
in
    getHexCorners

----------------------------------------------------

= Table.AddColumn(xy2latlon_2, "getHexCorners", each getHexCorners([latlon_L2], init_L2))        

Wait. Cant we do this in dax?

Sure you do! But unfortunately, the dax I used was a ton slower than Power Query. Now, I didnt pursue further to optimize this, but I think Power Query is a much better place to handle this than a calculated column, unless Brian Julius can use his super-ai stuff to make something that performs that much faster than the above power query. Is it a challenge Brian???

(oh and note, the dax below didnt have any init-style things baked in either. Not that it cant be done, but I just stopped trying further here)

WKT  = 
VAR _level = hextable[l]
VAR center_lat = 52.5055
VAR center_lon = 4.9590
VAR edge_length_lat = 174.38 / 1000 / 111.32/ sqrt(7)^(_level-9)
VAR edge_length_lon = 174.38 / 1000 / (111.32 * COS(RADIANS(center_lat))) / sqrt(7)^(_level-9)
VAR rotation = [angle]
RETURN
  IF(ISEVEN(_level),
    "POLYGON((" & 
        CONCATENATEX(
            GENERATESERIES(0, 6, 1),
            FORMAT([rotated_lon] + edge_length_lon * SIN(RADIANS([Value] * 60 +  rotation)), "0.000000") & " " &
            FORMAT([rotated_lat] + edge_length_lat * COS(RADIANS([Value] * 60 +  rotation)), "0.000000") & 
            IF([Value] < 6, ",", "")
        ) &
        "))",
        hextable[WKT])        

Closing word

This completes the code to find the hexagon corners on EP5! Once you have read it its pretty straightforward, but I wanted to give you the versions I came through in my journey. Perhaps it increased the learning experience a bit, perhaps you found it dauntingly boring and didnt even come to the end of this.

In any case, let me know if you like this better than the carousel or not or if you have some comments/remarks