Thursday, 25 January 2018

coordinate system - ArcGIS 10.0/10.1 - arcpy - automate projection tool with different datum or projection


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()

Answer



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...