Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
    
I left the Python and PostGIS tags off because I wasn't aiming for platform specific answers. I probably should've been explicit about that. –  canisrufus Jan 23 '12 at 22:25
    
Your solution's a pretty good one. Out of curiosity, what would the results of a Roeck test mean when applied to Alaska? That God used the Pacific plate to gerrymander the Aleutian Islands? :-) –  whuber Jan 23 '12 at 22:35
1  
@whuber You'd be surprised at the political clout Concerned Sea Otters for Fair Geological Formations have. ;) –  canisrufus Jan 24 '12 at 0:13
1  
Hi @whuber, now that I've changed jobs and will never, ever go back: When I was working for the Florida legislature, someone above me wanted to make the argument that the legislature can't draw compact districts because Florida is itself not compact. To that end, I was running compactness tests on all fifty states to "prove" it. To answer your question, "what would the results of a Roeck test mean when applied to Alaska": not a whole lot. But calling your boss mendacious and asking him not to waste your time is, for whatever reason, taboo, so I did it. –  canisrufus Jul 18 '12 at 20:10
add comment

1 Answer 1

up vote 1 down vote accepted

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()
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.