topic Re: Creating Lines of Graticule in Python in Python Questions

Creating Lines of Graticule in Python

EricShreve1 — Sat, 11 Dec 2021 12:12:57 GMT

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])

Re: Creating Lines of Graticule in Python

NeilAyres — Sun, 26 Mar 2017 12:11:31 GMT

Why not just use Create Fishnet to do this?

Create Fishnet—Help | ArcGIS Desktop

Re: Creating Lines of Graticule in Python

EricShreve1 — Sun, 26 Mar 2017 15:30:03 GMT

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

Re: Creating Lines of Graticule in Python

JoshuaBixby — Sun, 26 Mar 2017 15:50:48 GMT

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.

Re: Creating Lines of Graticule in Python

JoshuaBixby — Sat, 11 Dec 2021 12:13:00 GMT

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.

Re: Creating Lines of Graticule in Python

EricShreve1 — Sun, 26 Mar 2017 17:04:15 GMT

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.

Re: Creating Lines of Graticule in Python

JoshuaBixby — Sat, 11 Dec 2021 12:13:03 GMT

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.