Guess my EPSG

Let's say someone has a shapefile and they want to know the appropriate EPSG code for it. How do I do that??

So far I can see that most of the files have a descriptive string for example

"NAD_1983_StatePlane_California_II_FIPS_0402"

and they also have units

UNIT["Meter",1.0]

That's often all I need.

Here is the code that works for me...

#!/usr/bin/python
import sys,os
import fnmatch

lov = [
    ('"NAD_1983_StatePlane_California_II_FIPS_0402"',   'UNIT["Meter',  26942),
    ('"NAD_1983_StatePlane_California_II_FIPS_0402"',   'UNIT["Foot',    2226),

    ('NAD_1983_California_Albers', 'UNIT[Meter,', 3310),
    ('"NAD83 / California Albers"', '', 3310),

    ('"NAD_1927_StatePlane_California_I_FIPS_0401"', 'UNIT["Foot_US"', 26741),

    ('PROJCS["Custom",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",-4000000],PARAMETER["Central_Meridian",-120],PARAMETER["Standard_Parallel_1",34],PARAMETER["Standard_Parallel_2",40.5]', 'UNIT["Meter",1]', 3309),


    ('"Old_Hawaiian_StatePlane_Hawaii_3_FIPS_5103"', 'UNIT["Foot_US"', 3563),


    ('"NAD_1983_StatePlane_Nevada_West_FIPS_2703_Feet"', '', 102309),
    ('"NAD_1983_HARN_Lambert_Conformal_Conic"', 'UNIT["User_Defined_Unit",0.3048', 2913),

    ('"NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601_Feet_Intl"', '', 2913),
    ('"NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601"', 'UNIT["Foot"', 2913),
    ('"NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601_IntlFeet"', '', 2913),
    ('"NAD_1983_HARN_StatePlane_Oregon_South_FIPS_3602_Feet_Intl"', 'UNIT["Foot', 2914),
    ('"NAD_1983_HARN_StatePlane_Oregon_South_FIPS_3602_AltUnits"', 'UNIT["User_Defined_Unit",0.3047999902464]', 2914),

    # DOGAMI uses this and it's sort of a guess
    ('PROJCS["_MI_0",GEOGCS["NAD83",DATUM["D_NAD83",SPHEROID["Geodetic_Reference_System_of_1980",6378137,298.2572221009113]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",43],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_origin",41.75],PARAMETER["central_meridian",-120.5],PARAMETER["false_easting",1312335.958],PARAMETER["false_northing",0]', 'UNIT["Foot (International)",0.3048]]', 3644),

    ('"NAD_1983_Oregon_Statewide_Lambert"', 'UNIT["Meter', 2991),
    ('"NAD_1983_Oregon_Statewide_Lambert_Feet_Intl"', '', 2992),
    ('"OGICdefaultAG82"', 'UNIT["Foot"', 2992),

    ('"NAD_1983_HARN_Lambert_Conformal_Conic"', 'UNIT["Foot"', 2994),

    ('"NAD_1927_StatePlane_Oregon_North_FIPS_3601"', 'UNIT["Foot_US', 32026),
    ('"NAD_1927_Oregon_North"', 'UNIT["Foot_US', 32026), 

    ('"NAD_1983_StatePlane_Washington_South_FIPS_4602_Feet"',      'UNIT["Foot_US', 2286),

    ('"NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet"', 'UNIT["Foot_US', 2927),
    ('"NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602"',      'UNIT["Foot_US', 2927),

    ('"NAD_1983_StatePlane_Washington_North_FIPS_4601_Feet"', '', 2285),

    ('NAD_1983_UTM_Zone_11N', 'UNIT["Meter"', 26911),

    ('"WGS_1984_Web_Mercator_Auxiliary_Sphere"', 'UNIT["Meter"', 3857),
    ('"WGS_1984_Web_Mercator"', 'UNIT["Meter"', 3857),

    ('GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0]', 'UNIT["Degree",0.0174532925199433]]', 4152),

    ('GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]', '', 4269),

    ('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984"', 'UNIT["Degree"', 4326),

    ('GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]', '', 4326),

]

def find_epsg(prj):
    """ Point me at a PRJ string and I will attempt to tell you the EPSG code as an integer. """

    try:
        for tup in lov:
            (desc, units, epsg) = tup
            if prj.find(desc) != -1 and (len(units)==0 or prj.find(units) != -1):
                return epsg
    except:
        print tup
        exit(0)
            
    return None # huh?

#########################################################################

if __name__ == '__main__':

    matches = []
    for root, dirnames, filenames in os.walk('/g/GISData'):
        for filename in fnmatch.filter(filenames, '*.prj'):
            matches.append(os.path.join(root, filename))

    for prjfile in matches:
        fp = open(prjfile, 'r')
        line = fp.readline()
        fp.close()
        epsg = find_epsg(line)
        if epsg == None:
            print prjfile, line
            exit(1)
        else:
            print prjfile, epsg

    exit(0)