# Determine Unknown Coordinate Systems: Based on Spatial Extent Values

Discussion created by Playa on Feb 9, 2017
Latest reply on Apr 18, 2017 by curtvprice

I'm currently developing a Python Module to reproject shapefiles and CAD drawings (DWG, DXF) that have known and unknown coordinate systems. The code below deals with the known coordinate systems of the shapefiles of which there are 1500 shapefiles and 3600 CAD drawings. I realize that there are limitations in being able to identify the correct coordinate system of shapefiles\CAD drawings with unknown coordinate systems. I found the following two articles:

Projected Coordinate System Tables: Datum Arc 1960

Projected Coordinate System Tables: Datum WGS 84

I was hoping that by comparing the spatial extent of the shapefiles\CAD drawings with the extent of the coordinate systems (i.e. Arc_1960_UTM_Zone_37S, WGS_1984_UTM_Zone_37S) I would be able to determine the correct coordinate system. The first problem is that the extent is in decimal degrees and for some reason or other they are in WGS 84.

Projected Coordinate System Tables: values are based upon WGS84

If anyone has been able to come up with a workflow of identifying unknown coordinate systems based on the spatial extent of the shapefiles\feature classes\CAD drawings i'd truly appreciate any advice.

``'''Created on 03 Feb 2017@author: PeterW'''# import site-packages and modulesimport arcpyimport glob2import osfrom pathlib import Path# set argumentsinput_folder = Path(r"H:\Tanzania2")output_folder = r"H:\gdb"# set environment settingsarcpy.env.overwriteOutput = True# output coordinate system WGS_1984_UTM_Zone_37S \ 32737out_coord = arcpy.SpatialReference(32737)# function to create new geodatabasedef fgdb(input_folder, output_folder):    gdb = "{0}.gdb".format(input_folder.parts[-1])    arcpy.env.workspace = output_folder    arcpy.AddMessage("Creating {0}".format(gdb))    if arcpy.Exists(gdb):        arcpy.Delete_management(gdb)        arcpy.CreateFileGDB_management(output_folder, gdb)    else:        arcpy.CreateFileGDB_management(output_folder, gdb)    out_gdb = "{0}\\{1}".format(output_folder, gdb)    return out_gdbout_gdb = fgdb(input_folder, output_folder) # function to determine projection transformationsdef projections(shp):    spatial_ref = arcpy.Describe(shp).spatialReference    trans = ""    if spatial_ref.type != "Unknown" and "Cape" in spatial_ref.GCS.datumName:        trans = "Cape_To_WGS_1984_2"    if spatial_ref.type != "Unknown"  and "Arc" in spatial_ref.GCS.datumName:        trans = "Arc_1960_To_WGS_1984_3"    return spatial_ref, trans# function iterate over shapefiles, reproject to single coordinate systemdef reproject_shp(out_gdb, input_folder, out_coord):    # return shapefiles in Directory    arcpy.env.workspace = out_gdb    shp_pattern = "{0}\\**\\*{1}".format(input_folder, "*.shp")    for filename in glob2.iglob(shp_pattern):        try:            shp_records = int(arcpy.GetCount_management(filename).getOutput(0))            if shp_records > 0:                sr, transf = projections(filename)                if sr.type != "Unknown":                    fcs_name = arcpy.ValidateTableName(str(Path(filename).stem))                    out_dataset = "{0}\\{1}".format(out_gdb, fcs_name.title())                    arcpy.AddMessage("Projecting {0}".format(fcs_name))                    arcpy.Project_management(filename, out_dataset, out_coord, transform_method=transf)        except (RuntimeError, arcpy.ExecuteError) as e:            print(e)reproject_shp(out_gdb, input_folder, out_coord)``