Scripting to load only selected rasters based on an index shapefile

I realize I still haven’t gotten around to Part 2 of that other sequence yet. I will, but I hit a little nuisance over the weekend at work…well, it’s always been a nuisance, but a minor one. As we move to higher resolution imagery, pixel counts rise, file sizes or numbers increase, and I foresee that this is a nuisance that will grow, not go away, so I thought I’d see if I could make my life easier and my work more efficient–my time better spent, my employers’ money better invested–if I could figure out a way to load only the files I need to cover a certain area of interest (market area, administrative boundary,…?) without manually selecting them from a folder containing 6000 files. For whole rows, it’s not difficult…just select a range, but for some cases that’s inconvenient, like if you want a slightly irregular chunk down the eastern edge cut into six .ecws, stair-stepped along the eastern corner of the county…then you want sequences of uneven length from the right side of each row, also of uneven length, and you have to manually select these out of your window listing all the 6000 .tifs for the county to load into Global Mapper because reprojecting the whole county would tie up your machine for a length of time that you might like to reduce…

Where was I? Oh yes. It’s pretty easy to pop into ArcMap, overlay the two index grids, and Select By Location to figure out which .tifs I really need. I know from noodling around that it’s pretty easy to create a Search Cursor in Python to read all the filenames out of a field in my shapefile and spit them into a text file. This was a  n early babystep–I could at least then copy and paste 20 or so at a time from the text file into the Open>File dialog box, without having to manually look and remember or write down “10682 to 10473, 9986 to 9840, …” for each row. Just blithely cut and paste, 20 at a time. For 6 sids, at what, 300, 400 .tiffs at a time, maybe more? Yeah. Better, but not great.

So then the next step was figuring out how to automate the ingestion into Global Mapper of this information. Ideally, if ArcMap could talk to Global Mapper, I could have ArcMap tell Global Mapper which tiles I want to load.

Incidentally, I found that having Global Mapper import files, reproject, and export via script it works better and faster on my machines at work than it does if I have to wait for the display to refresh during the reprojection. And that was a direct offshoot of this.

I had noticed in Global Mapper’s File drop-down an option to “Run Script”, so I thought, hmm. If it understands Python, that would be handy, since I’ve already played with that a bit. It doesn’t. I looked into Global Mapper’s forum and Scripting Language Reference, and discovered in the process that the workspace (.gmw) is actually just a script that loads the files and sets the projection and display characteristics. It’s really pretty simple when you stare at it for a while. And exporting to a raster format is really pretty much just as easy as importing and setting up the workspace, except it says “export” instead of “import”.

So why not set up a tool in ArcMap to export a Global Mapper script (.gms), I thought. It appears the scripting in Global Mapper is such that I could probably set the whole thing up to run directly at a command from a tool integrated into ArcMap, but for now I figured I’d just run the scripts from IDLE and code paths directly in a text editor or shell, rather than deal with GUI inputs.

Anyhoo, so there you have it. I take two shapefiles and overlay them to find which features in the one layer intersect a certain feature or set of features from the other layer. Then I figured out what the Global Mapper script would need to say, which was basically:


followed by

IMPORT “C:\GIS\Blog\GeoTIFFs\” + file_name over and over.

This is what my Python script would write. Then if I wanted to, I could project it, say, from State Plane 6 to UTM 11 N, and export it from the script by adding a couple lines saying so after the import:


Projection UTM
Zunits NO
Zone 11
Xshift 0.000000
Yshift 0.000000



So, for simplicity and speed on my part, I would write a Python script that would take a shapefile as an input, create the header and import statements for the tiffs. To get this shapefile, I would do the selection in ArcMap and save the result. (Ultimately I could just create a layer and use the layer from a dialog box within ArcMap.) The shapefile I want it to read is cells on the tiff index grid selected because they intersect a selected grid cell on the coarser grid:

And the link between this shapefile and the Global Mapper script would be this Python script:

Then I can just open up Global Mapper, select Run Script, navigate to and select the .gms, and then hit go.


Evolution of a Python Script into a Script Tool: Part 1

The script I’ll be using in this 3-part series of posts is the latest iteration of the script in my last post. It goes back a few posts, and some may find it informative to trace it back to my effective, if not so user-friendly, initial version which required me to export a text file and so on. Some may find it informative; however, others may find it tedious. So let’s pretend you don’t feel like going back through those, and let’s start fresh. I will issue a little disclaimer–the script you see below did not emerge fully formed, but went through several steps before getting here, including bouncing back and forth through the tool validation and attempts to run the tool in ArcMap and running chunks of it in IDLE…

So here’s the script I said you’d see below, but I think maybe I should go back further. All the way back. I’ll do so after this little image:

If you’ve read the earlier posts, skip to the next paragraph. For those who haven’t, here’s the abstract: Where I work, we deal with imagery that is often tiled into a grid of tiffs, sids, or ecws. We have (or we create) an index shapefile that shows where each tile is in the grid. Sometimes I have wanted to rename the tiles for a variety of reasons. I figured that just as a field in the attribute table contains the current name of each tile, I could add a field and populate the names I want them to have. Then, using a Python script, I could read the values in these two fields for each record, and rename the files in a given folder from the existing name to the desired one.

So, the bare bones version of this script is in my last post. Basically, it takes four parameters: the index shapefile, the field name corresponding to the current file names, the field corresponding to the desired one, and the folder containing the files. Then the script creates a search cursor, and then goes through the records getting the values from each of the designated fields, and renames the files from the current name to the desired one.

Easy peasy, right? It does require the fields to have exactly what you want the names to be. Now, sometimes there’s a file extension present in the attribute table, sometimes not, so I would either have to customize the script by tacking on ‘+”.tif”‘ where the old and new names are, or simply add them into the attribute table with the Field Calculator. Not much trouble, but I thought if I make it a script tool, I could maybe give some options for dealing with those file extensions, and also with prefixes like “2011_UTM_LosAngeles_”, or like “M:\11_UTM_LA_tiffs\”, that sometimes precede the actual simple name of the tile. So I added a couple fields to my test tile index so I could play around with this a bit:

Now, before I really start fleshing out the bare bones, I want to make sure the bare bones work as a script tool, right? So I right-click in the ArcToolbox window, add my custom toolbox, and then right-click the toolbox and add a script. I give it a name, a label, and a short description, and I like to check “relative paths”. I hit “Next”, and navigate to my script. I hit “Next”, and then here’s where it can start to get a little complicated, depending on the script and the parameters. Keeping the script open in IDLE for reference, I fill in the parameters. I give them the following Display Names: “Index Shapefile”, “Field Containing Old Name”, “Field Containing New Name”, and “Workspace”. I set the data type for the first of these to “Feature Class”. For the next two, I choose “Field”, and for the last I choose “Workspace”. So far, so good. I highlight “Field Containing Old Name”, and in the lower Parameter Properties window, I scroll down to “Obtained from” and choose “Index Shapefile”. Furthermore, since I want to use text strings, I click next to “Filter”, choose “Field”, and then uncheck all but the “Text” box. I repeat this for “Field Containing New Name”, hit “Finish”, and now the script tool appears in my toolbox.

When I double-click, it opens up a dialog box just like any other ArcMap tool. I fill in the fields:

Hit go, and after a few seconds I can see in the ArcCatalog window that my files have been renamed! Yay!

In Part 2, I’ll define a function to modify the names, test it in IDLE, incorporate it into the script, add the necessary parameters, and add a simple little line to the tool validation tab.

Python script to rename image tiles (geotiffs, sids, ecws…) using the index shapefile

Elsewhere I posted a script that I wrote to allow me to do this same task, to rename imagery tiles based on attributes in the index shapefile. It wasn’t perfect, and it had a couple extra steps, but I wrote it pretty quickly and it worked well enough. If you want to, you can read more about it on that other post, or the example application linked from there.

Earlier today, however, I greatly simplified this script so that I can skip the step of exporting the attribute table to a text file and then go through all that parsing. Instead, this uses a search cursor. I wrote this one at home on a machine with ArcGIS 10. This current version basically runs from IDLE, but I’m going to adapt it for a script tool, which is why I have “arcpy.GetParameterAsText()” in there for the input parameters. To run it from IDLE, I just delete those and insert the shapefile, fields, and workspace into the code. After I get it running neatly from within Arc 10, I’ll probably make a 9.3 compatible version to run at work. Here it is:

Python script to rename files based on a shapefile attribute — example application

I thought it might be helpful to post not only scripts that I’ve created, but why and how I use them. As written, these scripts tend to be fairly short and simple, but also very specific to some task I do repeatedly at work, or that I suspect I may have to redo multiple times with slight changes, and which would be tedious to do manually. In short, I thought an example might provide a helpful context for the scripts themselves.

This example uses the script from the post titled “Python script to rename files based on a shapefile attribute”.

At work, I deal with aerial photography. This imagery may come to us in tiles, a grid of geotiffs, for example. Sometimes customers desire that we deliver the imagery to them cut into a different tiling scheme in order to keep consistency with their imagery from previous years, or with some other data layers or boundaries. Sometimes, when compressing an area to .ecw or .sid we may need to break up an area into some arbitrary tiling scheme to balance file size with image quality. In those cases, I find the script from my post “Python script to rename files based on a shapefile attribute” to be useful. I could also see this as possibly be useful with, say, a NED dataset from the USGS with a metadata shapefile.

Perhaps in another post, I’ll go into more detail about how to create a tile index for the imagery — for now I’ll just say in brief that it can be accomplished in a few ways — it can be done using ER Mapper (create a mosaic using the Mosaic Wizard, then create a Mosaic Tile Index Grid using the mosaic) , Global Mapper (select layers in the Control Center, right click, and choose Create Area Features from Selected Layer Bounds), and I think with a little creativity there is some way to do it ArcMap 10 as well.

Anyway…here are some steps and screenshots of how I might use this script:

First, let’s start with a raster dataset–I’ll use imagery –and an index shapefile of some sort.

The attribute table currently has a field called “TIFF” that contains the number of the tiff that fills the extent of the corresponding feature. Now, I want to replace this incomplete, sequential, bottom left to top right numbering scheme with one that is more convenient for our customers using rows and columns. I begin by adding three fields: “TILE”, “Row”, and “Column”.

I use the Select Features tool to select all features in a row, then use the Field Calculator to populate the “Row” field with “A” for the selected features. Rinse and repeat until all records have a value for row and column. Manually selecting rows and columns can be a little unwieldy for a large number of rows and/or columns — I may try to write something to automate this step eventually, but so far I haven’t found it necessary.

Once I have the rows and columns done, I use the Field Calculator to fill the “TILE” column. In this example I simply use the row and column separated by an underscore.

Now, if I change the label field to “TILE” instead of “TIFF” for the index shapefile, the new naming scheme is displayed.

Now I just have to get the names of the tiffs to match! First, from the Attribute Table for the index, I export to a text file. The way I have the script set up, this table has to be in the same folder as the tiles to be renamed.

Now I open the script in IDLE, and make sure the path, text file, and the fields corresponding to the old name and the desired new name are set.

Run the script, and the files are renamed to match the “TILE” field in the index shapefile.

At this point I’d probably load the newly renamed tiffs into ArcMap along with the index and make sure everything is good, then delete the “TIFF” field from the shapefile and LookupTable.txt from the folder to avoid confusing the customer. And, we’re done!

Some hints and resources

I realize I have not been terribly thorough in giving instructions with the scripts I have posted so far. They have been short and relatively simple, but in order to use them, I think you would have to be familiar with the process of running a Python script from IDLE and/or creating a script tool in ArcMap. I am using ArcGIS 10 at home, but at the office I have ArcGIS 9.3.1 installed, so the scripts and the methods by which I run them are different depending if I am at home or at work. I will try to remember to be specific about which I’m using, and give more details about how to run these scripts or import them to script tools and so forth in the future.

In the meantime, though, here are a few hints–

(1) If you see “import arcpy” near the top, I am using a computer with ArcGIS 10 and Python 2.6. If you see “import arcgisscripting” I am using a computer with ArcGIS 9.3.1 and Python 2.5.

(2) If you see lots of variables hard-coded, meaning things like “elevInput = ‘C:\\GIS\\TEST_FOLDER\\elev_10m'” then I have written the script to run from IDLE. If you see “elevInput = arcpy.GetParameterAsText(0)”, then I have written it to be used as a script tool.

(3) If you are saying, “huh?” or don’t know about scripting but want to learn, then I recommend checking out the following resources:

First, get familiar with ESRI’s online help if you aren’t already. For ArcGIS 10, go to and get familiar with the Geoprocessing Tool Reference (Professional Library > Geoprocessing > Geoprocessing Tool Reference). Also check out the bit on Creating Script Tools (and off of that, Customizing Script Tool Behavior) under Geoprocessing. For 9.3.1 Geoprocessing and Geoprocessing Tool Reference are both on the root at These are extremely helpful as they provide examples and details about the parameters.

Okay, now Python-wise, I recommend Mark Lutz’s book Learning Python ( Personally I bought the 3rd Edition because it was cheaper and covers up through 2.6, which is what Arc 10 uses.

Also, through Penn State you can access quite a lot of stuff from an online course for free at There are other resources worth a look at, while you’re at it.

And you might as well read through some of Dive Into Python:, and bookmark Okay, now don’t get overwhelmed– I think it was Lao Tzu who said what has generally been translated to “A journey of a thousand miles begins with a single step.” The journey may seem overwhelming, but you can surely handle a single step. And after that, take another when you are ready. And another. And we’re underway!

Python script tool to export each feature in a feature class to its own new feature class

Here is the script:

Selects each feature in a feature class and exports each to a new shapefile.
Created by John Speargas, 29 Jan 2012

# import modules
import arcpy
from arcpy import analysis, env

# input parameters
in_fc = arcpy.GetParameterAsText(0)
field = arcpy.GetParameterAsText(1)
env.workspace = arcpy.GetParameterAsText(2)

# create search cursor
rows = arcpy.SearchCursor(in_fc)

for row in rows:
where_clause = field + “='” + row.getValue(field) + “‘”
# name the output features based on the value found in the field selected above
analysis.Select(in_fc, row.getValue(field), where_clause)

This script is set up to be added to a toolbox in ArcGIS 10. I set up the script tool parameters as shown in the dialog box below. When you select the second parameter, the field based upon which the new features should be named, in the lower part of the dialog box click next to “Obtained from” and select the input feature class parameter.

Python script to rename files based on a shapefile attribute

The script is below. In another post, I show an example of how I use this script. It meets a specific need I encounter periodically at work. The application is pretty specific, but I think it is an example of how a simple script can augment GIS skills, so look for that post coming soon. Please note that if you copy this script from this page and paste it into IDLE you will probably have to go through and replace the quote marks (both single and double), as well as the two minus signs, as they seem to copy over incorrectly. This can be accomplished quickly by using “Edit > Replace”, copying the incorrect quotation mark into the “Find” field, and entering the correct one in the “Replace” field. Repeat for each incorrect character.

Okay, here it is:

Renames files in a target folder using a lookup table exported as a text
file from the attribute table of an ESRI shapefile.  The lookup table should
be in the same folder as the target files, or the code changed accordingly.
John Speargas, 2011

# Import module
import os

# Input folder containing files to be renamed
folder = “C:\\GIS\\TEST\\directory_to_be_renamed\\”

# Input lookup table
lookupTable = open(folder + “LookupTable.txt”)

# Read lookup table and close file
tableList = lookupTable.readlines()

# Input fields corresponding to the old and new names
old = “OLD”
new = “NEW”

# Read the lines in the table and strip away unnecessary characters
newTable = []
for line in tableList:
    newLine = line.strip()   
    newerLine = newLine.strip(“,”)  
    newestLine = newerLine.strip(“‘”)   
    almostFinalLine = newestLine.replace(‘”‘,””)  
    finalLine = almostFinalLine.split(“,”)

# Get position of old and new field names
header = newTable[0]
indexOld = header.index(old) – 1
indexNew = header.index(new) – 1

# Make lists corresponding the new and old values
oldList = []
newList = []
for item in newTable:
    if item != newTable[0]:
        valueOld = item[indexOld]
        valueNew = item[indexNew]
# For each file in the folder, parse the old name from (the path
# and-?) the extension
for file in os.listdir(folder):
    nameAndExt = os.path.splitext(file)
    oldName = nameAndExt[0]
    extension = nameAndExt[1]

#   Index the old name in the old list if it is present
    if oldName in oldList:
        indexOldName = oldList.index(oldName)

#       Locate the same index in the “new” list
        newName = newList[indexOldName]
#       Rename the file using the path, the new name, and the extension
        os.rename(folder + oldName + extension, folder + newName + extension)