diff --git a/qsomap.py b/qsomap.py index 1ec15f8..28a355a 100755 --- a/qsomap.py +++ b/qsomap.py @@ -103,6 +103,36 @@ def is_point_in_polygon(point, polygon): return True # odd number of intersects -> inside polygon +def polygon_centroid(polygon): + """ Calculate the centroid (center of mass) of the given 2D polygon. + + See https://en.wikipedia.org/wiki/Centroid#Of_a_polygon for the source + formulae. + + Args: + polygon: 2xN numpy array with x and y coordinates of N points. + + Returns: + (x, y): The coordinates of the center of mass. + + """ + # Separate the coordinates + x = polygon[0, :] + y = polygon[1, :] + + # First, calculate the area of the polygon + A = 0.5 * np.sum(x[:-1] * y[1:] - x[1:] * y[:-1]) + + # Calculate x and y of the centroid + C_x = np.sum((x[:-1] + x[1:]) + * (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * A) + + C_y = np.sum((y[:-1] + y[1:]) + * (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * A) + + return C_x, C_y + + def svg_make_inverse_country_path(doc, map_radius, polygon, **kwargs): # build a closed circle path covering the whole map commands = [f"M 0, {map_radius}", @@ -257,7 +287,8 @@ def svg_add_maidenhead_grid(doc, ref_lat, ref_lon, map_radius): lon = x * 2 * np.pi / N lat = np.linspace(-np.pi/2, np.pi/2, resolution) - x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, map_radius) + x, y = map_azimuthal_equidistant(lat, lon, + ref_lat, ref_lon, map_radius) points = np.array([x, y]).T + map_radius group.add(doc.polyline(points, **{'class': 'maidenhead_line'})) @@ -306,9 +337,10 @@ def svg_add_country_names(doc, simplegeodata, map_radius): w = np.max(x) - np.min(x) h = np.max(y) - np.min(y) - # FIXME: calculate center of mass - center_x = np.median(x) + map_radius - center_y = np.median(y) + map_radius + # align text at the center of mass + center_x, center_y = polygon_centroid(poly) + center_x += map_radius + center_y += map_radius kwargs = { 'class': 'country_label',