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?
Answer
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])
No comments:
Post a Comment