interactive map – test map up and running

Update, just in case anyone, anywhere ever reads this …

So, goings on — my friends Josh and Maxine were in town recently, and I got to go see their movie screening on Sunday at Santa Barbara, so that was very cool. I met a couple gents from Greenpeace, one of whom, at least, Josh and Maxine are have built a relationship with through the making of this film. I’ve pretty much decided that I’d rather work for an environmental nonprofit than just about anyone else, so I hope I’ll be in contact with those guys again. The movie is Musicwood, and it is nominated for a social justice award, so that’s exciting. Here’s the website: They’re headed to Missoula, Montana, Washington, D.C., and I think Seattle, Washington as well for film festivals in the not too distant future. Check it out.

And how about that interactive map? Well, I don’t think L.A. Waterkeeper would mind me saying that it is their website for which I’m making the map. I have a test map working, but I’ll wait to share until it’s better. But it’s going along pretty well. There were a few options, as mentioned in my previous post — OpenLayers, MapBox, ArcGIS Online — I decided to go with Google Maps API, because it’s pretty easy to use, has more tutorials and documentation than OpenLayers, seems to offer more in terms of adding links, icons, and so forth compared to ArcGIS Online, and there shouldn’t be any licensing issues for a nonprofit. I didn’t really look into MapBox as much — it does look pretty good, but really Google Maps just stood out to me as probably the best blend of being easy to learn, relatively quick to get set up, and though their basemap must be used, it can be styled using the Style Wizard to look very different, very uncluttered. I can add layers as KMLs, which allows for placemarks and line styles and all those sorts of things to be established within the code of the KML pretty easily. It’s a little tricky to muddle into, because the various GIS programs I’ve tried all do different things when exporting to KML — ArcMap allows you to set the colors, but it dumps all the attributes, and transparency must be added after export. QGIS did kind of a weird thing — the file I used had these simple but nice tables that showed in balloons when you click on a feature in Google Earth, and it kept all the attributes, but then in the Google Map none of that happens. It’s good to know the info is being retained within the KML, but I would have to restructure it to get it to display the way I want. Global Mapper so far is the best I’ve tried — it doesn’t have all the lovely qualitative color ramps of ArcMap, but it does keep the color, the transparency, and the attributes. It gives the features placemarks, which pop up a balloon containing a rather ugly and rudimentary little table, but I’m thinking maybe I can borrow the table style tidbits from the QGIS export as a template to make those tidier. So far, so good.

I think it’s going well, and I think the folks at Waterkeeper seem pretty pleased with what little I got going so far, and I think we’ll have it up on their website in a month or so, after which it will be all the sort of tweaks and refinements to improve it. Anyways, I wish I had more time, but my actual paying job has been keeping me busy. Plus there’s that pesky thesis I’m supposed to be working on.

interactive web map

I’m trying to create an interactive web map as part of an internship. As Yoda once said, “Do or do not, there is no try”. I know I can produce something that will be a useful display of information involving maps. One way would be to create a bunch of maps at predetermined extents and give the user the capability to navigate through a larger reference map to maps that are progressively zoomed in. I could do this by tiling, or by just making maps according to my own arbitrary choices of what looks good and conveys information. The bonus is that they have nothing right now, so anything I produce is an improvement. At the very least, it should be possible to turn layers off and on. One means of doing this is to add even more maps to the catalog above, the same extents but with different combinations of layers turned on and off. Now, this is starting to be a lot of maps to create, which means updating or changing the information would be a bit of a nuisance, but if I keep a map document handy I could just run through the bookmarks with a certain layer combination and export. I might even be able to script that. I should keep world or prj files with those, or just always use the same projection so I know what it is without one, so I could add those jpgs or whatever to a GIS. Ideally the layers would be queryable, and with permission, editable, but I consider that a goal to achieve later; first I need to see what’s involved in making it simply interactive on a web page. I’m thinking ArcGIS Explorer Online and Google Maps API will give me the capability to do the most interactivity with the most user-friendly interface, provided I’m okay with their logo (I am) and their basemap options underneath. I’m also looking at things like MapBox and TileMill, and OpenGeo which stacks PostGIS, GeoServer and Open Layers along with a couple others. I think the last one might be the best, as it actually involves a suite containing the database, the GIS, the server, all together, while the former I believe allows you basically to make maps and serve them, but I guess the data just needs to be in shapefile, KML, GeoTIFF, csv, or some other common format. So maybe I could just go back into a project and replace a “site” layer with an updated one, and how I store it and update could be Arc or Global Mapper or QGIS or OpenJump or whatever. Well, it’s Christmas eve, and I’ve got things to do. Cheers,


Lead Paint

Okay, what the hell? I was trying to keep this blog more or less professional, but I think it’s a pretty dry read so far, and I look up there and see my name in the address bar, and I think, well, that doesn’t seem right. So, I’ll add a little personal touch now, I think.

Recently, my friend Josh ( introduced me to this guy Adam (, and we played some music. If you scroll down there on Adam’s page you can see the MUSIC heading, and if you click Lead Paint (2008-Present), the following page will open: In the address bar, add “071” to get You can hear what Adam mixed from the recording of that session.

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.