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.

I am getting stuck with a task that is slightly different from the posts I found so far on this same website.

Assumption: all layers have a defined coordinate system

Problem: I would like to automate a python script (using arcGIS 10.1 but 10.0 is fine as well) that projects a list of vector layers (some of which with different coordinate systems) to a common coordinate system, given by a target vector layer. The main issue I would need help with is how to automate something this:

1) check if the coordinate system of the current vector layer is the same as the target layer

2) if not, check if the datum is the same

3) if it's the same datum, then project the current layer to the target layer without any datum transformation

4) if, instead, they do not have the same datum, pick a pre-defined transformation and include that in the projection function

Any help would be much appreciated! thanks!

Here is my sample code:

import arcpy
from arcpy import env

arcpy.gp.overwriteOutput = True
env.workspace = 'C:\\Temp'
Target_lyr = 'water_Broward.shp'

try:
 desc = arcpy.Describe(Target_lyr)
 TargetSR = desc.spatialReference

 for inputFC in arcpy.ListFeatureClasses('*'):

  if inputFC != Target_lyr:
   inputFC_SR = arcpy.Describe(inputFC).spatialReference

  if inputFC_SR.Name == TargetSR.Name:
   continue
  else:
   #check if current layer and target layer are the same type
   #thus if they are both 'projected' or 'geographic' 
   if inputFC_SR.type == TargetSR.type:

    #HELP STARTING FROM HERE...
    #Check if both have same datum..

    #if they ARE, then
    #outFC = inputFC + '_prj.shp'
    #arcpy.Project_management(inputFC, outFC, TargetCS)

    #if datum is NOT the same use
    #transformation = for example "WGS_1972_To_WGS_1984_1"
    #arcpy.Project_management(inputFC, outFC, TargetCS,transformation)

   #if current layer and target layer are NOT same kind 
   else:
    #HELP
    #repeat as above same steps as above to check the datum 

except:
 print arcpy.GetMessages()
share|improve this question
    
Nearly identical question asked just a few hours ago: gis.stackexchange.com/questions/49086/… –  blah238 Jan 25 '13 at 23:17

1 Answer 1

up vote 4 down vote accepted

some relevant functionality below. all of which is new at 10.1.

with projected coordinate system, you have access to the GCS object which is where the datum info lives. if you're working with a GCS object already, just access the relevant datum properties (datumName, datumCode) directly.

from_sr = arcpy.SpatialReference('NAD 1983 HARN UTM Zone 11N')
to_sr = arcpy.SpatialReference('NAD 1927 StatePlane California VI FIPS 0406')

print("from sr datum: " + from_sr.GCS.datumName)
print("  to sr datum: " + to_sr.GCS.datumName)

# list transformations valid for San Diego county
outlist = arcpy.ListTransformations(from_sr, to_sr,
                                    arcpy.Extent(444450.2212, 3599832.1877,
                                                 585727.9387, 3707930.3429))

print("valid transformation: " + outlist[0])

The code above will give you these outputs

from sr datum: D_North_American_1983_HARN
  to sr datum: D_North_American_1927
valid transformation: NAD_1983_To_HARN_CA_S + NAD_1927_To_NAD_1983_NADCON

The ListTransformation's extent is optional but you could pass in you data's extent for best transformation. This is the same functionality used in gp's Project tool (to pick default transformation) as well as map's data frame. Documentation link

share|improve this answer
    
Thanks a lot gotchula! So you mean that in 10.0 I would have to find a different workaround? Any ideas? –  user20038 Jan 25 '13 at 23:45
    
at 10.0 it's more involved yeah. i'd look at parsing the output of sr.exportToString() to get the datum/gcs. as far as finding the transformations, there is no way to do that at 10.0 as far as i'm aware. –  gotchula Jan 25 '13 at 23:52
    
gotchula I tried your suggestion, but I still have a problem...if I get the extent of my target layer (with: desc = arcpy.Describe(Target_lyr).extent and the .XMax etc.) and then use that inside the arcpy.ListTransformations definition...in position outlist[0] it gives me one that is good for Alabama...which does not make sense since my target layer is defined over Broward county, FLorida...any idea why? Thanks! –  user20038 Jan 28 '13 at 16:12
    
your Target_lyr is the in_sr (as per the documentation)? beyond that i've no idea without looking at your data. If you want to bundle it with code i can have it looked at it. –  gotchula Jan 29 '13 at 1:05
    
Hi, my Target_lyr is the variable to_sr based on your first answer...in the arcpy.ListTransformations() I use the arcpy.Extent taken from Target_lyr...at least that is what I understood I had to do based on your answer..how can I bundle data with code? do you mind sending me an email to ftonini@ufl.edu so that we can do all that via email? Thanks! –  user20038 Jan 29 '13 at 16:19

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.