I am sure I am re-inventing the wheel here but I couldn’t find any ready-made scripts to do this. Plus, I thought, it was high time I get my feet wet with Python. So- behold of my first python script(!). Its quite simple- it basically takes an input linear shapefile and converts it to a new shapefile with M values. If a field containing the measured line length is given, it will use it to calibrate the M values – otherwise it will use the map (geometry) length so the M value of the last vertex of each line will equal the length of the line.
This is the first of a bunch of ‘useful’ Avenue scripts that I used quite often in the past and converted to Python. Will eventually convert the others and learn more about Python and arcpy in the process.
The script is listed below for your viewing pleasure. It is heavily commented so it should be self explanatory. Or, you can simply download it from here
import arcpy, math, os # Input feature class # inFeatures="C:/mapdata/osm_greek_data/projected/road_sample.shp" # Fieldname containing the measured length for each line # if null, then the map length will be used mlField="meas_len" # Output feature class settings out_path = "C:/mapdata/osm_greek_data/projected" out_name = "road_sampleM.shp" geometry_type = "POLYLINE" template = inFeatures spatial_reference=inFeatures has_m = "ENABLED" has_z = "DISABLED" # Create the new shapefile arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference) # # Create a list of fields to use when copying field values from input to output # fields = arcpy.ListFields(inFeatures) fieldList = [] for field in fields: fieldList.append(field.name) # Loop through input features rows=arcpy.SearchCursor(inFeatures) # Output features outFeatures=out_path+"/"+out_name # Open insertcursor # outRows = arcpy.InsertCursor(outFeatures) # Create point and array objects # pntObj = arcpy.Point() arrayObj = arcpy.Array() shapeName=arcpy.Describe(inFeatures).shapeFieldName #Get shape field name featType=arcpy.Describe(inFeatures).shapeType # Get feature type isMeasured=arcpy.Describe(inFeatures).hasM # Check if its a polylineM feature already # Record counter i=1 # Loop through input shapefile for row in rows: print "Processing row %d " % i feat=row.getValue(shapeName) noparts=feat.partCount partno = 0 currentM=0 # Loop through geometry parts for part in feat: pntCounter=0; segmentDistance=0 prevX=0 prevY=0 # Loop through points in part for pnt in feat.getPart(partno): if pntCounter == 0: # First point prevX=pnt.X prevY=pnt.Y pntObj.X=pnt.X pntObj.Y=pnt.Y pntObj.M=0 else: curX=pnt.X curY=pnt.Y # Calculate segment length segmentDistance = math.sqrt(math.pow((curX-prevX),2) + math.pow((curY-prevY),2)) # Use the measured distance if mlField <> "": segmentDistance=segmentDistance/feat.length*row.getValue(mlField) # Add previous M to segment length currentM=currentM + segmentDistance pntObj.X=pnt.X pntObj.Y=pnt.Y pntObj.M=currentM prevX=pnt.X prevY=pnt.Y # Add measured points to array arrayObj.add(pntObj) pntCounter=pntCounter+1 # Create new row # outFeat=outRows.newRow() # Copy attributes to new row # for fieldName in fieldList: if fieldName not in ("FID", "Shape"): # ignore FID and Shape fields outFeat.setValue(fieldName, row.getValue(fieldName)) # Assign array of PointMs to new feature outFeat.shape=arrayObj # Insert feature outRows.insertRow(outFeat) # Clear array of points # arrayObj.removeAll() i=i+1 # Remove cursor from memory del outRows del rows
