Python – Minimum Bounding Circle of Geometry Crossing the 180th Meridian

geometrypostgispython

Background: I've been given a task where I have to calculate the minimum bounding circle of Alaska (for the purpose of running a Roeck test). I'm using the Census' state level geometries available here.

The problem: Alaska's minimum bounding circle gets calculated "incorrectly". There is a string of islands located between ~172 and 180 degrees longitude, while the rest of Alaska is from -180 to -129 degrees longitude. This results in my minimum bounding circle covering the whole globe. While it is technically a correct minimum bounding circle of the given points on a euclidean plane, it clearly doesn't represent the true minimum bounding circle of Alaska. (I'm using PostGIS's ST_MinimumBoundingCircle.)

The solution: The way I made it work was to

  1. convert Alaska to well known text
  2. subtract 360 from all of the longitudes between 170 and 180 degrees
  3. Put Alaska back into the database.

(The code will be in an answer.)

This is a dissatisfying solution because it's not generalized (it won't work in all cases), and it seems harder than necessary.

My question: What are the options for dealing with this problem? Has someone written a more intelligent version of the MBC tool that will find the "true" MBC instead of the over-large one? Is there a simpler way to transpose the geometries than transposing individual points?

Best Answer

This is how I solved the problem:

import psycopg2
import re

def replace(matchObj):
    # Group 0 is all groups, group 1 is the first match contained in
    # parentheses. Group 2 is the whole float value.
    value = float(matchObj.group(2))
    value = value - 360
    return matchObj.group(1) + str(value) + " "

def main():
    # connect to postgresql
    conn = psycopg2.connect("conn info")
    cur = conn.cursor()
    cur.execute("select st_astext(the_geom) from tl_2010_us_state10 \
                WHERE name10='Alaska'")

    # Get query results
    rows = cur.fetchall()
    alaskaText = rows[0][0]

    # The regular expression finds a positive float between 170 and 179.9, 
    # inclusive. Matches any character not a negative sign to ensure that the
    # longitude value is positive, then it matches the rest of the number.  
    # To ensure that the number matched is longitude, we check that it is
    # followed by a space (the format being "longitude latitude,".
    matchText = re.compile('([^-])(1[67][0-9]\.[0-9]*) ')

    # Substitute longitude values
    alaskaText = matchText.sub(replace, alaskaText)

    # insert into dabase and commit transaction
    cur.execute("insert into tl_2010_us_state10 (the_geom, name10) \
                VALUES (st_geomfromtext('%s', 4269), 'Alaska2')" %(alaskaText))
    conn.commit()
    conn.close()

if __name__ == "__main__":
    main()
Related Question