0

I'm hoping to create a script that creates a new shapefile for each record in a table. My python experience is limited, and it has been a while since I used it, so I'm having a hard time getting started. My table has fields for ID, longitude, latitude, distance and bearing.

For each record, I need to:

  • create new shapefile with a filename that matches the ID field
  • create a series of points within the newly created shapefile at 1 meter increments along a line that is defined by the longitude/latitude (this is the start point), distance and bearing fields.

Here is my initial code:

import arcpy
table = "C:/folder/project.gdb/transect_table"
fields = ["TRANSECT_ID", "LATITUDE", "LONGITUDE", "Heading", "Distance"]
shp_path = "C:/folder/output_shapefiles"
with arcpy.da.SearchCursor(table, fields) as cursor:
    for row in cursor:
        arcpy.CreateFeatureclass_management(shp_path, cursor[0])

This bit of code is creating the collection of empty shapefiles with filenames being pulled from the ID field of the table.

Can anybody suggest an approach to populating the shapefiles with points?

Based on the input from radouxju, I've updated the code to the following (note that I'm creating feature classes in a gdb as opposed to shapefile; my shapefiles were getting locked when creating the points but feature classes seem to avoid this issue):

import arcpy
import math

path = "C:/folder/Output_FeatureClass.gdb"
table = "C:/folder/Transects.gdb/table"
fields = ["TRANSECT_ID", "LATITUDE", "LONGITUDE", "Heading", "Distance"]
with arcpy.da.SearchCursor(table, fields) as cursor:
    for row in cursor:
        arcpy.CreateFeatureclass_management(path, cursor[0], "Point","","","", "C:/folder/WGS1984.prj")
        with arcpy.da.InsertCursor(path + "/"+  cursor[0], ["SHAPE@XY"]) as insertcursor:
            for i in range(row[4]):
                xy = [(row[2]*i*math.sin(row[3]),row[1]*i*math.cos(row[3])),]
                insertcursor.insertRow(xy)

This is working pretty well, as all of my feature classes are getting created and are populated with the appropriate number of points. However, the distance between the points is not the desired 1 meter (it's closer to 15,000 km!). I suspect that I might have an issue with mixing units, as the input positions are in degrees, and the offset between points should be in meters, but I'm not quite sure.

Any suggestions?

  • Here's some starting direction: You'll need a loop to go through all of the features, and I suggest two functions: One function which takes in the start-point, bearing and distance and returns a collection (the series of points based on the equation of the line.) the second function can write output based on the loop-index, the selected table record (for the name) and the point collection returned from the first function. – JasonInVegas Jan 13 '15 at 1:39
  • 1
    I think this is a multi part question where the first part should be to figure out how you are going to read values from the table. If the table is in a geodatabase or if it is in dBase format then I would say to use a Search Cursor (preferably from arcpy.da - but the availability of that will depend on your version). Perhaps figure out how to read from a table in this question and then look at writing features (using an update cursor) separately. – PolyGeo Jan 13 '15 at 1:40
  • 4
    How far have you got? i.e. what have you tried already.. As PolyGeo said there's at least two questions here; have a read of resources.arcgis.com/en/help/main/10.1/index.html#//… for the iteration and if you're still unsure of where to go next post a specific question. Of course if you have some working code please edit your question and include some so we can tell where you're up to, otherwise we may be covering topics you're already familiar with and just wasting words. – Michael Stimson Jan 13 '15 at 1:55
  • Your question has been placed On Hold because, as mentioned in the Tour we are looking for one question per question. To get it re-opened I recommend using the edit button to focus it down to how to use ArcPy to read (and perhaps print values) from your table. Once you have your code to do that then it will be easier to research/ask how to write data to a new shapefile separately. – PolyGeo Jan 13 '15 at 4:46
  • The next bit "create new shapefile with a filename that matches the ID field" should be easy - just use arcpy.CreateFeatureClass_management() - if you can get that included and working in your code before anyone offers an answer then you may yet be able to start looking at how to write your point features within this question. – PolyGeo Jan 13 '15 at 7:01
1

I have not tested but it should be something like below. I assume your heading is 0 toward North and in radians, and that you are in a projected coordinate system with the units in meters (ideally Mercator).

import arcpy
import math
table = "C:/folder/project.gdb/transect_table"
fields = ["TRANSECT_ID", "LATITUDE", "LONGITUDE", "Heading", "Distance"]
shp_path = "C:/folder/output_shapefiles"
with arcpy.da.SearchCursor(table, fields) as cursor:
    for row in cursor:
        arcpy.CreateFeatureclass_management(shp_path, cursor[0], "POINT")
        with arcpy.da.InsertCursor(shp_path + "/"+  cursor[0], ["SHAPE@XY"]) as insertcursor:
            for i in range(row[4]):
                insertcursor.insertRow((row[2]+i*math.sin(row[3]),row[1]+i*math.cos(row[3],))

if your heading is in degree, you'll need a conversion from degree to radian (math.radians(row[3])

  • Thanks radouxju. I tried your code but I'm getting an invalid syntax error related to the SHAPE@XY (at least that is where my cursor jumps to). Any thoughts? – bryanhi Jan 13 '15 at 22:16
  • " " forgotten, sorry – radouxju Jan 13 '15 at 22:19
  • Next problem...I'm getting an error that 'sin' is not defined. I've tried a few basic things like importing math module at start of script, and typing calcs as math.sin and math.cos, but still no success. – bryanhi Jan 13 '15 at 23:08
  • Actually I think its correct to import the math module, and add the math.sin and math.cos. The other thing that was missing from your script was that a trailing comma on the tupple containing the new coordinates (see: support.esri.com/en/knowledgebase/techarticles/detail/40604). However, something seems to be wrong with the math. Instead of 1 meter increments, I am getting increments of like 50,000 kilometers. I actually cant even get arcmap to measure properly since they are so far apart. Any thoughts on fixing the math? Thanks for the help. – bryanhi Jan 14 '15 at 2:20
  • thank you for testing. There was indeed something wrong with the math (multiplying by increment instead of adding the increment). the code has been updated. – radouxju Jan 14 '15 at 6:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.