Creating Lines of Graticule in Python

1897
6
03-25-2017 10:30 PM
EricShreve1
New Contributor II

I am stuck with what step I need to take to create a loop statement to automate the process of creating lines of graticule in python/arcpy. What the end result is to create a shapefile that draws a polyline every 10 degrees for lines of latitude/longitude. An example would be the first line drawing at (-90,180), (90,180) and second line drawing at (-90,170), (90,170) until it gets to -180. Here's my current status of my code. 

# Import ArcPy
import arcpy

# Create empty shapefile
arcpy.CreateFeatureclass_management('C:/Python_Testing/Week5/Week5/Data', 'line1.shp','POLYLINE',spatial_reference=arcpy.SpatialReference(4326))
# Create cursor to insert data
cursor = arcpy.da.InsertCursor('C:/Python_Testing/Week5/Week5/Data/line1.shp', ["SHAPE@"])

# Writing a multipart feature from an array of points from the same cursor


firstPart = arcpy.Array([arcpy.Point(180, -90),
                         arcpy.Point(-180, -90)])
secondPart = arcpy.Array([arcpy.Point(180, -90),
                          arcpy.Point(180, 90)])

array = arcpy.Array([firstPart, secondPart])

polyline = arcpy.Polyline(array)
cursor.insertRow([polyline])
Tags (2)
0 Kudos
6 Replies
NeilAyres
MVP Alum

Why not just use Create Fishnet to do this?

Create Fishnet—Help | ArcGIS Desktop 

EricShreve1
New Contributor II

Unfortunately, the class I am taking requires submiting a python document.

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

Since this is for a class, I will demonstrate creating the longitude lines and you can modify the code to generate latitude lines:

SR = arcpy.SpatialReference(4326)

long_lines = []
for x in xrange(-180,190,10):
    line = ((x, -90 + y/10.0) for y in xrange(180*10 + 1))
    long_lines.append(
        arcpy.Polyline(
            arcpy.Array(arcpy.Point(*coords) for coords in line),
            SR
        )
    )

arcpy.CopyFeatures_management(long_lines, "in_memory/long_lines")

The code above creates longitude lines every 10 degrees.  Instead of densifying afterwards, the code creates lines that are densified at 0.1 degree increments.

EricShreve1
New Contributor II

Joshua, 

That's exactly what I need! Formatting the "For Loop" was confusing me. Like you mentioned in the previous response I would have gone the route of using AGOL to query the data but for scripting purposes, I needed to make the code. 

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

No problem Eric.  As much as grabbing a standard Esri data set from AGOL makes sense in most cases, there is also value in learning to gin up the data yourself.

For Line 8 above, or what looks like Line 7.5 in my browser, I pass a generator expression to ArcPy Array to build the polyline.  If you are learning Python, I strongly encourage you to read up on comprehensions (list, dictionary, set) and generator expressions.  Comprehensions are typically some of the most performant ways to generate lists, dicts, etc...., and I think most people find them very readable as long as the logic doesn't get too complex in them.

As an example, my original code can be re-written using a list comprehension:

SR = arcpy.SpatialReference(4326)

long_lines = [
    arcpy.Polyline(
        arcpy.Array(arcpy.Point(x, -90 + y/10.) for y in xrange(180*10 + 1)),
        SR
    )
    for x in xrange(-180, 190, 10)
]

arcpy.CopyFeatures_management(long_lines, "in_memory/long_lines")

Whether the list comprehension is more readable than my original code is arguable.  For me, I find working and reading comprehensions very natural, but I have been working with the language for years.  I think for people new to Python, they take some getting used to.

JoshuaBixby
MVP Esteemed Contributor

Is there a reason you have to generate your own?  You could just download World Latitude and Longitude Grids from ArcGIS Online and use a definition query to select the 10 degree increments you are interested in using.

Say for whatever reason you want or need to create your own and from scratch using ArcPy.  One of the biggest problems you will run into with your current approach is that your lines will not project properly in most world projections.  There are a few world projections, like Web Mercator, where the lines you create will line up nicely with actual latitude and longitude, but for most that won't be the case because you haven't created a curve but a straight line with only two points.

There are several ways to address the projection issue.  One is to create densified lines instead of a line with only two points.  You can either densify the line when constructing it or make the lines and then run Geodetic Densify on the final feature class.

Another approach would be to build true curved lines.  This is a bit tricky because you would have to use Esri JSON to define the lines since ArcPy Geometry constructors don't work with curves.  You would also need to make some decisions on whether to use simple circular arcs as approximations or elliptical arcs.  If the latter, you would need to make assumptions about major and minor axis.  Although I would try the true curves for kicks, creating the lines and densifying is much more straightforward.

Regarding existing geoprocessing tools, I have found Create Fishnet to be clunky when working with hemispheric, and especially global, grids.  There is also the Make Grids and Graticules tool, but that has a bit of a learning curve.

The data team at Esri has already makes some of these choices and done the work for folks, I encourage you to think about downloading and using the data I referenced above if possible.