Convert a Polyline shapefile to PolylineM using arcpy

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
Advertisements

Shp2ora and Ora2shp: Utilities for importing and exporting shapefiles to Oracle

Yes, another one of those. Its not exactly a novel idea as there are a few converters around for moving shapefile data to and from Oracle:

  • If you are an ESRI/ArcSDE user you have the sde2shp and shp2sde commands.
  • Oracle provides two options: a command line program, shp2sdo to load shapefiles into Oracle, and its java equivalent through Mapbuilder
  • Most (all?) GIS’ that can read Oracle Spatial data, include options to save the Oracle layer to a shapefile.

Plenty of choice you probably think. Well, I was still not happy. First of all I didn’t always have ArcSDE on all machines I was working on, secondly installing Mapbuilder means you need to go through the 18.2MB download, plus it will only LOAD shapefiles into Oracle and not the other way around and thirdly, I didn’t want to go through the whole process of downloading a fully-blown desktop GIS (say QGIS) just for the sake of converting a spatial table into a shapefile. Not to mention that you have to configure the client GIS first in order to connect to Oracle BEFORE you can even display the layer – so you can then export it. Too much hassle for a relative easy task.

So I decided to create my own little command-line utilities to do just that. You can download the full source code and executables from the Box.Net widget on the left side of this post (ora2shp.rar).

After compiling the project, you will find two executables under the Debug folders of Ora2shp and shape2ora projects: Ora2shp.exe and shp2ora.exe respectively,

Ora2shp syntax is: ora2shp <username/password>@dbalias> <spatial_table_name> <PK_col> <shape_col> <shapefile> ["optional_where_clause"]

If you use the where clause you should NOT put the WHERE keyword.

Shp2orashp2ora <username/password>@dbalias> <spatial_table_name> <shape_col> <shapefile> <srid>

If you don’t specify a SRID it will default to null. If the spatial table already exists records will get appended. If it doesn’t it will get created and will also create the oracle metadata records record (USER_SDO_GEOM_METADATA table) and spatial index.

It should deal with all shapes (multipart, Z, M) APART from multipatch.

Note that you will need Oracle client 10.2 installed to compile the programs as-is. Otherwise change the reference to the Oracle client of your choice. You will know if you have the wrong Oracle client version if you get an error about Oracle.DataAccess not being the same version as Oracle Client.

When running the programs any errors should appear on the console and the full stack trace will be written to orashp_err.log file, located at the same folder as the executables.

Note that  performance is not great- especially when creating the shapefile. It took something less than 15mins to create a point shapefile of around 66K records. (around 5M of shapefiles). But feel free to improve it and I would very much appreciate to get back to me if you do!