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 have a large data set X(x1-xn) from which I must select consecutively the attributes A(a1-an) from X.fieldsSRC or X.fieldDST and project the selections.

I have a functioning arcpy script that does the following:

<inside loop incrementing n>
#proj_sr is defined here but omitted
where_clause = fieldSRC = n OR fieldDST = n #simplified for webparser

#create a new feature class based on selection
arcpy.Select_analysis('input', 'temp1', where_clause)
#project selection
arcpy.Project_management('temp1', 'temp2', proj_sr)

#PERFORM OPERATION ON PROJECTED SELECTION

#append selection, and project to destination SR
arcpy.Append_management('temp2', 'final_output')

#we must delete rows from source files otherwise we will treat them each twice
#create layer file in order to select and delete specific rows
arcpy.MakeFeatureLayer_management('input', 'layerfile', where_clause)
#delete rows selected above in layer file
arcpy.DeleteFeatures_management('layerfile')
<end of loop>

This works. In my mind there should be a manner to do the following which would be more efficient.

<inside loop incrementing n>
#proj_sr is defined here but omitted
where_clause = fieldSRC = n OR fieldDST = n #simplified for webparser

#create layer file in order to select and delete specific rows
arcpy.MakeFeatureLayer_management('input', 'layerfile', where_clause)
#project ACTUAL selection rather than a complete shapefile as was done earlier
arcpy.Project_management('layerfile', 'temp2', proj_sr)

#PERFORM OPERATION ON PROJECTED SELECTION

#delete rows selected above in layer file
arcpy.DeleteFeatures_management('layerfile')
<end of loop>

Obviously the arcpy.Project_management() function does not operate in the second code example. It's a rather large data set ('input') and recreating the feature class based on the selection each time does not seem efficient/necessary. It works but it is not fast.

Does anyone have a suggestion?

share|improve this question
    
I think you should pursue your second code example further. MakeFeatureLayer accepts an SQL expression and Project accepts a feature layer so I am thinking it should work unless there was an code error there when you tested that. Perhaps get just those few lines of code out into a separate test and report back on what issue is arising from it. –  PolyGeo May 10 '12 at 3:56
    
The problem I'm encountering using the desired method (the second code block) is that Project is selecting the full source data set rather than the layer selection. The arcpy.MakeFeatureLayer_management operation has created the layer 'layerfile' from 'input' successfully and selected the desired x rows. When I pass 'layerfile' to arcpy.Project_management it projects all the rows from 'input'. –  Cyrille May 16 '12 at 7:49

3 Answers 3

up vote 1 down vote accepted

UPDATED ANSWER

The code below should work using the file geodatabase that Cyrille posted to the Esri Discussion Forum.

import arcpy

prjString = "PROJCS['World_Mercator',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],UNIT['Meter',1.0]];-20037700 -30198300 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision"

arcpy.CopyFeatures_management(r"C:/Support/Cyrille/MFL_Test.gdb/FlightPaths",r"C:/Support/Cyrille/FlightPaths.shp","#","0","0","0")
arcpy.MakeFeatureLayer_management(r"C:/Support/Cyrille/FlightPaths.shp","YVR_layer",'"src" = ' + "'YVR'")

arcpy.Project_management("YVR_layer", r"C:/Support/Cyrille/MFL_Test.gdb/YVR_Proj", prjString)

What I found is that my original code made the layer on a shapefile while Cyrille was making it on a file geodatabase feature class. To workaround what I think is either a MakeFeatureLayer bug or a Project bug what I did was to copy the feature class to a shapefile before performing those two steps on it. That way the where clause is respected.

I suspect the workaround will not be desirable for Cyrille so the next step will be to see if this bug has been reported previously and to report it if not.

ORIGINAL ANSWER

I just ran the test below which uses MakeFeatureLayer_management to apply an SQL expression to the underlying dataset (of 31 polylines), and then uses Project_management to project that layer (of one polyline). The results is a feature class of one projected polyline. So it looks like what you are trying to do in your second code example should work, unless your where_clause is not evaluating correctly. Consequently, my recommendation is to double-check that your where_clause is working.

import arcpy

arcpy.MakeFeatureLayer_management("C:/Temp/test.shp","test_Layer",'"ID" = 555')

arcpy.Project_management("test_Layer","C:/Temp/test.gdb/projectedFC","PROJCS['GDA_1994_Australia_Albers',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',132.0],PARAMETER['Standard_Parallel_1',-18.0],PARAMETER['Standard_Parallel_2',-36.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]","#","GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
share|improve this answer
    
I created a super simple example using my data and posted the geodatabase and code on the esri forum. I agree and believe it should work. It clearly is not though. I used the python window so I can clearly see that the 'where_clause' is functioning correctly. Give it a spin. –  Cyrille May 18 '12 at 14:51

Have you tried using a SearchCursor and setting the spatial reference on it so it reprojects the results on the fly? Code below straight outta the help:

import arcpy

desc = arcpy.Describe("d:/base/data.gdb/latlongbnd")

# Create search cursor. Use the spatial reference object from the
#   described feature class so geometries are returned in decimal degrees.
rows = arcpy.SearchCursor("d:/base/data.gdb/buildings", "", desc.spatialReference)
share|improve this answer
    
I do not need to iterate through a list. I need to select 1-200 rows of spatial data from a file containing ~40,000 elements and project just the selection. While iterating may not be slower than a selection (although I hope it is), the projection of all the elements may well be. I would have to do some timing tests but I am unsure that this would be faster. –  Cyrille May 7 '12 at 8:21

SearchCursor takes a "where" clause - it's the best way to deal with a subset of records.

share|improve this answer
1  
Could you please explain the answer with detail? –  Sunil Apr 6 '13 at 11:46

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.