8

I am using ArcGIS 10.4.1 and Python 2.7.10.

I have a mxd file of a rather large area which contain hundred of layers. I wish to work in a specific sub-area in that workspace, where only a few layers will be relevant.

I was looking for a way to output a list of the layers which are spatially overlapping with the current view (sub-area) in ArcMap. The layers could be visible or not, as checked (or not) in the workspace. The mxd is up to date with this view.

Is such a thing possible?

I am using as a basis the following code, which output all the layers in the workspace:

import arcpy

mxd = arcpy.mapping.MapDocument(r"C:\path\to\file.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
print arcpy.mapping.ListLayers(mxd,"", df)

EDIT: a practical example to illustrate would resemble to this: Let's say I have dozen of layers from topographical maps and ground features (water, river, topoline-level, roads, electric line, settlements, glaciers, reservations, highways, airports, and forests) at a country scale.

I zoom in a very remote forested area, where only the layers forest, river, topoline-level are overlapping the view. All other are not relevant. I am looking if it is possible to arcpy those forest, river and topoline layers in a list and discard the other.

4
  • 1
    OK What if you had the situation of a donut, all your data is around the edge but none of it is visible in the current view but the extent of the dataset covers it? As @Aaron says update your question with more information rather than responding in the comments.
    – Hornbydd
    Jan 26, 2017 at 17:11
  • @Hornbydd and Aaron thanks for feedback - I clarify based on your respective comments in the original question post
    – marsisalie
    Jan 26, 2017 at 17:18
  • So in the donut scenario (e.g. settlements around your forest but not in the view) you would not want to report this layer?
    – Hornbydd
    Jan 26, 2017 at 17:22
  • @Hornbydd this is correct ! Since they are around and not in the view those would not be required in this scenario
    – marsisalie
    Jan 26, 2017 at 17:33

2 Answers 2

9

A DataFrame object has an extent property, the extent property/object can be used in the basic spatial relationship methods of contains, within, equals, overlaps, touches and disjoint as well as distanceTo, see here for more details on those methods.

With that in mind you can loop over the layers you retrieve from calling arcpy.mapping.ListLayers(mxd,"", df) and retrieve each layers extent by calling the getExtent method on each layer. The spatial relationships discussed above can be tested between the extent.polygon object of each layer and the dataframe to determine if they are indeed within the extent of the dataframe. If you are using data driven pages the extent of each page can be retrieved with an mxd object with this call mxd.dataDrivePages.dataFrame.extent.

It is important to note that the extent of a layer often has a larger footprint that than the actual vector within the layer (e.g. the polygon, points, raster, etc.), so depending on the spatial relationship you choose to test it may include more layers than are actually visible. Here is an example (thanks @KHibma) of how you would use a spatial relationship method of the DataFrame with a Layer object:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for l in arcpy.mapping.ListLayers(df):
    if not df.extent.disjoint(l.getExtent()): 
        print l.name
7
  • if you dont want that sample code, feel free to remove
    – KHibma
    Jan 26, 2017 at 17:40
  • @KHibma much appreciated. Perhaps I should have included it.
    – GeoSharp
    Jan 26, 2017 at 17:43
  • I was working on the code as an answer, but you spent more effort explaining the processes than I would have :)
    – KHibma
    Jan 26, 2017 at 17:46
  • 1
    You're probably better off going with 'not disjoint' or 'overlaps'. 'Touches' will almost never be true.
    – phloem
    Jan 26, 2017 at 17:57
  • 1
    People, seriously you are awesome - this is working perfectly. I had to remove/fix unresolved layers (layers in error due to wrong path) and thus with the disjoint directive it does the work. In addition, I learned about the Dataframe extent property which will be useful for upcoming projects. Warm thanks.
    – marsisalie
    Jan 26, 2017 at 19:21
7

As a quick and dirty method, you could compare extents, as already answered. However (also mentioned), that doesn't necessarily indicate whether a feature is visible in the current view (e.g. if you zoom in far enough within a feature class' extent, at some point you will not see any features).

There's a definite performance hit if you need to find visible features. You can get at it using several methods, but in some way you need to perform an overlay operation at the feature level. You could use: Spatial Join (as below), Intersect, Select Layer By Location, or search cursor through all features checking 'overlaps' against the df.extent, among others.

mxd = arcpy.mapping.MapDocument("CURRENT") # map
df = arcpy.mapping.ListDataFrames(mxd)[0] # data frame
arcpy.env.addOutputsToMap = False # disable add outputs
for lyr in arcpy.mapping.ListLayers(df): # loop through layers
    try:
        sp = arcpy.SpatialJoin_analysis(df.extent.polygon,lyr,r'in_memory\sp',"JOIN_ONE_TO_ONE","KEEP_COMMON",'',"INTERSECT") # perform overlay
        if int(arcpy.GetCount_management(sp).getOutput(0)) > 0: # check if overlay found any features
            print '{}: visible'.format(lyr.name)
        else: # no features found
            print '{}: not visible'.format(lyr.name)
    except:
        print '{}: not a feature layer or some other problem'.format(lyr.name) # not feature layer, e.g. WMS
3
  • Will definitely consider along other solutions
    – marsisalie
    Jan 26, 2017 at 19:23
  • 1
    Agreed, from a "completeness" factor, this is better. But the question/task is one that could be solved with "quick and dirty" or "slow and more correct". I like both answers. In the end its a matter of solving said problem in whatever way you're happy with. For that reason they both get +1
    – KHibma
    Jan 26, 2017 at 20:18
  • Sure, all fair points. However, I assume the main reason why someone would do this is to avoid drawing unnecessary, large datasets (e.g. countrywide base data for a site map), which are all the more likely to intersect a data frame if only the extent is considered.
    – phloem
    Jan 26, 2017 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.