Time Out
Dear readers,
This blog is on hiatus because my job situation has changed. I have taken a position as a programmer. At the moment I’m developing a database-driven web application with a combination of VB.NET and Javascript. I’m not working directly with GIS right now, although that is going to come back in to the mix eventually. I’m excited about this career move because I have been wanting an opportunity to try my hand at coding for quite some time. I would like to keep blogging, but the subject and the target audience is going to be different, and I need to decide how I am going to handle that.
If you’re still using ArcGIS version 9.3 or previous you might find something of use in my old entries. If you’d like to follow the adventures of an entry-level programmer then stay tuned!
Till later,
Aubrey
Editing OpenStreetMap data with QGIS and Merkaartor
Problem: You want to contribute to OpenStreetMap from your desktop instead of from the Potlatch web tool.
In my last entry I discussed the process of downloading OpenStreetMap data using Quantum GIS (QGIS). This is how the downloaded data looked in the vicinity of the North Plateau Loop in Monte Sano State Park:
In this entry, I will cover two different methods I used to update this data to include my GPS track of the entire loop. I initially attempted to do this with QGIS.
Method #1: QGIS
When you first activate the QGIS OpenStreetMap plugin, the editing widget should appear on the right side of the screen. If it doesn’t, turn it on by going to View –> Panels –> OSM Feature. You should use the widget instead of the standard toolbar to identify and edit OSM features. If you click the OSM identify tool and hover over features, they will turn red and you can see their tags.
Shown above (in red) is one of the footway bits that fell near the North Plateau Loop, along with my GPS tracks (in magenta). I decided to keep the two segments that were already there, along with their foot=permissive and highway=footway tags. I added a name=North Plateau Loop tag to both of them by typing in the <new tag here> area.
Then I pressed the Create Line button.
The editor will automatically snap to existing features. I began at one the end of the segment by the parking lot and traced along the length of my track, stopping to incorporate the second segment along the way. Then I added the three tags to the new lines I had created.
Once I was happy with my edits, I tried to upload them to the OpenStreetMap server. At this point in the process, if you do not yet have an OSM account, go to https://www.openstreetmap.org/user/new and create one. Then, go back to QGIS and press the Upload OSM Data button.
The next screen is a summary of your changes. Add an optional comment to explain what you did, then enter your OSM account username and password, and press the Upload button.
If all goes well, you will get a message saying your upload has completed successfully. I, instead, got a python error message “OperationalError: SQL logic error or missing database”. According to the Quantum GIS forum, this has happened to some other people, but it doesn’t happen to everyone. You may get lucky. If not, you can try another tool such as Merkaartor.
Method #2: Merkaartor
While QGIS is a full GIS software package, Merkaartor has been developed solely for the purpose of working with OpenStreetMap. I found the interface less than intuitive.
After opening the program, my next step would typically be to add a basemap that moves me to the right geographic location. However, there is no Add Data button like I would expect. Instead, there is a list of layers already in the table of contents that have no source. You need to right-click on the top one to point it to some data. You can choose from a pre-loaded list of TMS (tile map service) or WMS (web map service) layers, or add a raster in GeoTIFF or “Walking Paper” (JPG, PNG or BMP) format.
I had trouble getting Merkaartor to load my georeferenced trail map, so I went the map service route. Once you pick a layer (I liked OSM Mapnik), you can use the mouse wheel to zoom and the right mouse button to pan to your area of interest. Then, press the green arrow Download button to get the current OSM layers in that area.
In the next screen that comes up, you would choose to download the Current View, which is selected by default. Alternatively, you could choose to download “From the map below” and do the same zooming and panning to your area of interest within the smaller window.
Merkaartor’s way of displaying the Monte Sano OSM layers is a bit different. It shows every vertex (node).
(click to see full size)
When you hover over a group of nodes, their color will change to a pink highlight. You must click on them to select them. They will then turn royal blue and their tags will be displayed in the Properties window on the bottom left of the screen. In order to select the line instead of the nodes that make it up, you need to move the mouse until the pink highlight turns thicker. After you click, the resulting royal blue selection symbol will also be thicker. Then you will see the tags shown above. You can edit them by typing inside the fields, just as with QGIS.
Go to File –> Import to load your GPS track .gpx file. It will be symbolized in a way that is difficult to distinguish from the OSM data. There is no way to change the symbols. I recommend turning off every layer that you don’t need to make things easier to see. (Do this by clicking on the eye symbol next to the layer.) When you are ready to begin editing, go to the Create menu and select Road.
I started, once again, at the bit of footway near the parking lot.
A couple cool things will happen at this point that didn’t happen in QGIS. First, the program will snap your new line not only to the ends of the existing segments, but also to the vertices of your GPS track. This makes it easier to follow it exactly. Second, the program will automatically carry all of the tags associated with the existing segment over to your new road. You will not have to go back and add them later.
When you have finished digitizing, press the ESC key to exit creation mode. At this point, if you want, you can polish up your work by splitting and merging segments. Simply hold down the CTRL key to select multiple segments, then choose the function to apply from the Road menu.
I found the ability to do that very useful because there was a bit of North Plateau Loop that didn’t belong in the original segments. Also, I worked in pieces so I could zoom in closer, and I liked being able to merge everything into one at the end.
When you are ready to upload your changes, you will need to go to the Tools –> Preferences –> Data tab to enter your OpenStreetMap username and password. Then press the megaphone Upload button. Review and summarize your edits in the next screen.
If you’re me, this time it will work! There is now a perfect North Plateau Loop on OpenStreetMap.org for all to see.
(click to see full size)
My parting thoughts are thus:
I wish I had been able to complete the upload in QGIS, because I find the program more pleasant to work with, and I like the idea of using a single GIS package to do everything. However, Merkaartor does offer some advanced OSM editing functions that make for a cleaner final product.
In summary, these are the main benefits of each tool:
QGIS
- Full-service GIS package
- Allows for changing of symbol representation
- Allows for raster georeferencing and display
- Better interface for those familiar with ArcGIS
Merkaartor
- Allows for merging and splitting road segments
- Snaps to GPS trace vertices
- Copies tags from attached segment into extension
- Uploads OSM edits without issues
Tags: free, gps, gpx, OpenStreetMap, Web Map Service
Working with GPS tracks in Quantum GIS (QGIS)
I am beginning to explore the world of free GIS software. QGIS has been recommended to me as one of the best out there. As an ArcGIS user, I found it easy to learn because the features are very similar. I am also impressed with the number of things it can do. I am going to write about just a small sample of them:
- Georeferencing raster imagery
- Loading OpenStreetMap data
- Viewing GPS tracks and waypoints (.gpx files)
My goal when I started this project was to contribute some new information to OpenStreetMap. Since I arrived in Huntsville, I have had a lot of fun walking around Monte Sano State Park. I also recently obtained a GPS unit that I’ve been using for car navigation, and I thought I could use it to record the path of one of the hiking trails.
On my first visit to the park, the gate attendant gave me a trail map that looked close to this.
(Click to see full size)
A blue line surrounding the parking area, called the North Plateau Loop, looked to me like a good place to start. It passes some lovely scenic viewpoints.
I wanted to check if this trail has been added to OpenStreetMap yet. I decided to georeference the trail map and then download the OpenStreetMap data that falls within its bounds. Georefrencing within QGIS is a little different than what I am used to because it uses a separate window.
To start the process, select Georefrencer from the Plugins menu. Then, press the Open Raster button in the new Georefrencer window, and navigate to your image. Mine was montesano.jpg.
Press the Add Point button to start adding control points. When you click on the image, a box will come up that allows you to enter x,y coordinates or load the coordinates from the map canvas (your other window).
I loaded some vector street lines into my other window that I got from the City of Huntsville GIS Department’s FTP site. I georeferenced to the street intersections. When you’ve added enough control points–at least three, but more is better, go to the Settings menu to open the Transformation Settings.
The QGIS User Guide explained that Thin Plate Spline is a good transformation type to use when you have a lower-quality original, because it will introduce local deformations. Cubic is a good resampling method for smoothing things like scanned maps where we don’t need to worry about maintaining precise grid cell values. The only output format available is GeoTIFF. Chose the save location, then press ok. To initiate the georeferencing process, press the green start arrow.
If you checked “Load in QGIS when done”, the image will show up in your map canvas table of contents, and you can check your work. Everything will stay loaded in the Georefrencer window the way you left it, so it’s easy to add or remove control points and try again.
Now, it is time to load OpenStreetMap data on top of my map. The OpenStreepMap plugin is included with the core QGIS software, but it isn’t turned on by default. Go to the Plugins menu and select Manage Plugins to check it on. This will add the OSM buttons to the toolbar and will open OSM widget to the right of the map window. (I will talk more about the widget in my next entry.)
I zoomed so the trail map filled the screen, then pressed the Download OSM data button.
QGIS will grab the current extent of your map window. If the area is too big, it will tell you and you will need to zoom in. Select a place to save your .osm file and then press the Download button.
QGIS will automatically symbolize the layers similar to how they look on the OpenStreetMap website. When I looked at the data, I could tell that someone had added a couple footways near the bottom half of the North Plateau Loop, but they were not connected or complete.
So, I headed back to Monte Sano to walk the trail with my GPS device. I have a Garmin Nuvi 255w. The OSM wiki has instructions for how to collect traces with this device, but I didn’t find that I needed to do all that. My Nuvi has a Trip Log feature that is actually always recordings paths. It saves them to a standard .gpx file. You don’t need to tell it to start doing that, but you can turn on the lines so you can see it working. To do that, press Tools on the opening screen, then Settings, and then choose Show under Trip Log.
Now, whenever you move, you’ll see a cyan line being built behind you. If you want to start anew without all the clutter from your previous trips, press Tools, then My Data, then Clear Trip Log. I told my Nuvi to navigate me to the front gate of Monte Sano. Then I turned off my car, but left the Nuvi on and took it with me. I walked the trail and watched the line grow.
When I got home, I plugged the Nuvi into my computer and grabbed Current.gpx from the Garmin/GPX folder.
I added it to QGIS along with the OSM layers I had downloaded previously. To add a GPS track, press the Add vector Layer button. You will need to select the GPX file type from the long list in the browse data box.
Then, a window will come up asking you to choose which GPX sublayers to load. Tracks is the only one I need for this project. Waypoints holds your saved favorite places, and track_points holds the time when each point along the track was recorded.
I changed the symbol of my track line to a thicker bright magenta so it would show up more clearly. You can get to the symbology settings by right-clicking on the layer and selecting Properties, just like you would in ArcMap.
You can see that the GPS recorded my driving, parking, and walking movements all in the same layer. I think that the wiki instructions I linked to above might have allowed me to isolate just the trail, but then I would have had to do an extra file conversion step, and it is pretty easy for me to tell which is which.
My plan next was to edit the OpenStreetMap data that I downloaded to include this trail, and then upload my changes to the OSM server. I will compare the process of doing that with QGIS and Merkaartor in my next entry.
Tags: free, georeference, gps, gpx, OpenStreetMap, trip log
Hot spot analyis of water bodies in Madison County
Problem: You would like to determine where spatial features have unusually high or low concentrations
At the beginning of October, I will be moving to Huntsville, Alabama aka “Rocket City.” I have been thinking a lot about the city that is going to be my new home lately, and I have also been learning a lot about spatial statistics.
The city of Huntsville falls mostly within the county of Madison. The Census Bureau provides county-wide water body layers, as well as the area of water in each census block.
I thought I could use these data layers to help me decide where I should look for a house if I want to live near a lot of water. I am going to do Hot Spot Analysis on both of the layers.
Note: Hot Spot analysis is better suited for use with cultural, biologic or economic data, rather than physical landscape features. But I am working with what GIS files I could find, and I thought it would be interesting nonetheless.
Before proceeding, I need to project my layers out of geographic coordinates. I chose
NAD 1983 StatePlane Alabama East because that is used by the City of Huntsville GIS Department.
Let’s start with the water bodies. I am going to get a count of how many water bodies fall within each census block. To do this, first I converted the layer to points through ArcToolbox’s Feature to Point tool. Then, I right-clicked on the census blocks layer to do a spatial join with the water body point locations. This function will add a Count_ field to the attribute table of the output.
Next I ran hot spot analysis on that Count_ field. In 9.3, this is in ArcToolbox under Spatial Statistics Tools –> Mapping Clusters.
In the result, the red “hot” zones show where lots of different water bodies are clustered close together. The beige zones are average and the blue “cool” zones show where the water bodies are dispersed especially far away.
Now I will run the same analysis on the area of water within each census block. If you look at the attribute table, you’ll see that we’re given the land area and the water area. Those numbers aren’t very useful the way they are. I would rather have a percentage of the block area that is water. So I added another field: TotArea, that is the sum of the land area and water area. Then I added a Ratio field, that is the water area divided by the total area.
Hot spot analysis of this ratio field produces a map that shows where blocks with a high percentage of area covered by water are clustered close together.
Not surprisingly, the hottest areas are near the Tennessee River. A bit of a different story than the first map, which highlighted areas with lots of little lakes.
In these first runs, I let ESRI take the defaults for my parameters. I can improve my results if I pick the search distance on my own. The Spatial Statistics Tools –> Analyzing Patterns toolset has two tools that I used for this purpose.
The first: Spatial Autocorreclation (Morans I) will determine the correct scale for analysis by selecting the distance at which the factors promoting clustering are most pronounced. The tool takes a distance value as input and returns a z score as output. The absolute value of this z score indicates the amount of clustering. We want to find the distance that produces the z score farthest away from zero (either positive or negative).
I ran the tool a bunch of times, using distances from 0.25 – 1 mile in increments of 1/4 mile.
I found that 0.75 miles (3960 feet) produced the largest z score for the water area layer.
Note: A z score between -1.96 and +1.96 indicates a random distribution of features, which can be explained by the null hypothesis. Hot spot analysis is useless in this case, because there is no pattern to analyze. Any score outside that range indicates an unusually clustered or dispersed distribution of features that is influenced by some factor (or factors) other than random chance.
I knew to stop incrementing the distance when my z scores peaked and then started going back down. When I tried the Spatial Autocorrelation tool on the water count layer, this didn’t happen. The z score just kept going up. So I switched to the Multi-Distance Spatial Cluster Analysis (Ripleys K) tool.
I tried a bunch different settings on this tool, and got the clearest results with the ones below. Note I used the water body centerpoints as my input and employed a boundary correction method to reduce errors at the edge of the study area.
(click to see full size)
The tool produces a Diff K value. Like the z score, it indicates the amount of clustering, and we are looking for the largest number.
My best distance would be 7920 feet (1.5 miles) where Diff = 1763.49. Using these distances in the “Distance Band or Threshold Distance” box in the Hot Spot analysis tool produces maps that are more statistically meaningful.
(click to see full size)
This map is a merge of the hot and cold spots I generated for both the water body count and percent water area, using the best threshold distances in each case. It is telling me to move downtown if I want to avoid mosquitoes and to the outskirts if I want natural beauty.
Tags: cluster, hot spot, spatial autocorrelation, statistics
Creating a list of keywords from the attribute table: New & Improved Version
Problem: You need a better way to summarize the contents of a dataset without wasting time typing a long list.
Back in March, I posted about creating a list of keywords from the attribute table. I knew that:
The solution I have come up with is not perfect. I keep thinking there must be a better way.
I had explained how to make a list of the administrative area types in the Natural Earth Admin 1 shapefile, using ArcMap’s “Summarize” tool. The end result was a text file that had extra quotation marks, commas, and other characters to remove.
There is a better way — with Python! I now have the power to structure the text file the way I want it from the start. Here’s the entire code:
import arcgisscripting
gp = arcgisscripting.create(9.3)
report = file(r"C:\gis\AdminTypes.txt", "w")
shape = r"C:\gis\10m-admin-1-states-provinces-shp.shp"
cur = gp.SearchCursor(shape)
row = cur.Next()
adminList = []
while row <> None:
adminList.append(row.ENGTYPE_1)
row = cur.Next()
adminList = list(set(adminList))
adminList.sort()
for admin in adminList:
report.write(admin + "\n")
report.close()
print "Finished"
Change only the portions marked in red to make this work with any data layer. They represent, in order from top to bottom:
- The output text file
- The input data layer
- The name of the attribute table field to summarize
The text file generated by this code is a list of each unique entry in the ENGTYPE_1 field (all duplicates removed) — one per line, in alphabetical order. It’s the exact same thing you would see in the table of contents if you symbolized by category on that field. Only now it’s in a form you can copy and paste into other places, like a metadata record.
Python makes everything better.
Tags: attribute table, keyword, list, python, report, text file
Comparing attribute tables with Python
- Problem #1: You are merging layers that have different attribute table structures and you don’t want to lose any important information.
- Problem #2: You want to use Python to return a list of shapefiles that match a search phrase.
- Problem #3: You want to use Python to compare two lists and return the differences.
The following will solve all three!
I learned a lot from taking Introduction to Geoprocessing Scripts Using Python. When I returned to work, I got the grand idea that I was ready to write a script that would make one looming task more manageable.
As is always the case, the solution seemed easier in my head than it turned out to be when I sat down to write it. The ESRI course had given me just enough information to be dangerous. I had most of the pieces I needed to put together this puzzle, but not all. Those missing pieces are what I would like to share with you today.
The Background
I had created a geodatabase feature class and I needed to append and bunch of shapefiles to it. All of these shapefiles had attribute tables that were close to each other, but not exactly alike. Some had extra fields. I needed to make a list of all of those extra fields so that I could add them to the feature class before I started loading shapefiles.
The Solution
The Python arcgisscripting module has a ListFields() function that forms the basis of my solution. But before I could get to that point, I needed to make a list of the shapefiles whose fields I wanted made into a list.
Depending on your work flow, that might be all the shapefiles in a particular directory. For me, it was all the shapefiles in a particualr directory that matched a file name search phrase. One example: I wanted to merge together all the shapefiles that had “bridge” in the name. My teacher introduced us to the os.walk() function, but I needed to modify it slightly to do exactly what I wanted. Here’s how:
# import modules and create the geoprocessor object
# (some of these you won't use until later)
import arcgisscripting, os, fnmatch, sys
gp = arcgisscripting.create(9.3)
# returns a list of shapefiles in the specified folder
# that match the specified search pattern
folder = r"F:\StagingArea\shapefiles"
pattern = sys.argv[1]
shpList = []
for path, dirs, files in os.walk(folder):
for filename in fnmatch.filter(files, pattern):
shpList.append(os.path.join(path, filename))In my setup, the search folder is hard-coded and the search pattern is a variable that the user inputs. This worked for my purposes because I was searching through the same folder on a lot of different file names. Modify according to your purposes. Note: The search pattern needs to be in the form *bridge*.shp. You can use * wildcards and you need to put the .shp on the end.
The code so far has only created the list. To make it visible, use this:
gp.AddMessage("Shapefiles returned:")
for shp in shpList:
gp.AddMessage(shp)Both segments of code assume you will be running the script from within ArcToolbox. Use a print statement if you will be running it from within IDLE or PythonWin. Note that this list returns not just the shapefile names, but their full paths. This makes the list usable as input into the next function.
The code to create a list of attribute field names is pretty simple. Here’s how it looks for my single target feature class (the one I’m going to be appending everything to):
targetLayer = sys.argv[2]
targetList = gp.ListFields(targetLayer)
targetNames = []
for field in targetList:
targetNames.append(field.Name)To create a list of attribute field names for a list of shapefiles, you add another loop. Here’s how it looks using the shpList I created earlier.
x = 0
fieldNames = []
for shps in shpList:
fieldList = gp.ListFields(shpList[x])
for field in fieldList:
fieldNames.append(field.Name)
x = x + 1Now we have two lists. The first list is all the field names in the target layer and the second list is all the field names in every layer we want to add on to the target. What we really need is the difference between these two lists. In other words: we need to know which fields appear in the second list and not in the first list.
That sounds pretty complicated, but guess what, you can do it with two lines of code:
diffList = [fname for fname in fieldNames if not fname in targetNames]
diffList = list(set(diffList))
The first line compares the two lists. The second line’s list(set()) function removes any duplicates and alphabetizes the final list of differences.
Finally, to display your results:
gp.AddMessage("Missing fields:")
for field in diffList:
gp.AddMessage(field)My finished code does a few extra things, like:
- eliminating unmatching geometry types from the shapefile list
- eliminating empty fields from the field name list
- listing the field length and type in the final output
Perhaps I will talk about those in a later entry. What I’ve shown here are the parts that were the most difficult to figure out.
Tags: attribute table, geodatabase, list, python, search, shapefile
Using the command line to convert between metadata formats
Problem: You have a lengthy metadata record that you want to preserve, but it has been created using a different stylesheet than your standard.
Just a few days before I was going to take a course in python scripting, I noticed the command line button on the ESRI toolbars.
It is there in ArcMap and ArcCatalog. I wonder why I’d never noticed it before. Turns out you can do a lot of cool things with it. The portion of every Help document labeled “Command line syntax” is talking about things you can type into the command window. Plus, I learned in my class later that with ArcGIS 10, it will be renamed to the Python Window because we will be able to run python scripts directly from it! But I am getting ahead of myself.
There are certain things you can only do from the command window, and today I will talk about one of them: converting metadata from one format to another. The two prominent metadata standards are FGDC: Federal Geographic Data Committee and ISO: International Organization for Standardization 19115.
If you want to use data from an organization that follows a different standard than you do, you might need to do a conversion. For example, the metadata for World Wildlife Fund’s Global Lakes and Wetlands Database is in FGDC format. In order to view it in ArcCatalog, I need to have an FGDC stylesheet selected.
With the right stylesheet, all of the pertinent information is displayed in the metadata window.
If I switch to an ISO stylesheet, what I see is not nearly as helpful.
You could consign yourself to switching between stylesheets, but I suspect many people work in places that insist all of the metadata for their collection of geospatial records be in the same format. If you need to do a conversion, there is a faster way than typing all of the information back into the metadata creation wizard.
ESRI has provided a set of four metadata translators with the ArcGIS desktop install. They should be located in the \Program Files\ArcGIS\Metadata\Translator folder. You can run these translators from the command line window. First, pick the one you need.
- FGDC to ISO 19139 (FGDC2ISO19139.xml)
- ESRI-ISO to ISO 19139 (ESRI_ISO2ISO19139.xml)
- FGDC to ESRI-ISO (FGDC2ESRI_ISO.xml)
- ISO 19139 to ESRI-ISO (ISO19139_2ESRI_ISO.xml)
Then press the command line button, and type ESRITranslator into the window. If you press the spacebar, a tool tip will pop up that tells you the arguments you need to enter next. In this case it will say
ESRITranslator <source> <translator> {output} {logfile}
Input within carrots <> is required; input within brackets {} is optional. For me, <source> is the path to the Global Lakes and Wetlands metadata .xml file. <translator> is the path to the ESRI Translator I want to use ( FGDC to ESRI-ISO ). And {output} is the path to where I want to store the new converted .xml file. I’m not going to specify a logfile.
My code looks like this:
(click on image to see full size)
Notice a couple things here. My cursor is in the path to the translator file, and <translator> is bolded in the tool tip. This is intentional. If I move the cursor to other arguments, the bolding will change accordingly. Use this to make sure your syntax is correct.
Also notice that I have surrounded the path to the translator file in quotes. Quotes are optional, but you must use them if you have spaces in your path name. Otherwise the space will be misinterpreted as the beginning of the next argument. Going through and checking what is bolded where will help to catch this sort of thing.
Press enter to run the code. Then, press the Import Metadata button and browse to your new .xml file.
Now it looks how ISO should!
Tags: command line, fgdc, iso, metadata, xml
Doing network analysis
This is post #3 in a three part series. It is meant to be read after Creating a Network Dataset.
Now that we’ve done all the hard work building our network dataset, it is time to have some fun putting it to work. Most presentations would cover this part first to lure you in. I might have done things backwards. If I haven’t scared you away yet, and you’ve brought along a working network dataset, let’s explore what Network Analyst can do with it.
I will give examples of the first three Network Analyst tools: Route, Service Area and Closest Facility. Each of these revolve around facility locations paired with the network dataset. The City of Portland GIS Catalog doesn’t have a Points of Interest or Buildings layer, so I used ESRI Streetmap USA data for the point locations.
Route Tool
The Route solver finds the best path between any number of stops. Let’s pretend we are tourists who want to see a bunch of sites downtown and need some help planning our day so we’ll spend the least amount of time driving. I opened up ArcMap and added my new network dataset and the ESRI Recreational Areas point layer. Zooming to the Portland area shows 12 sites on both sides of the Willamette River.
Check on the Network Analyst toolbar, and then select New Route from the drop down menu. Also, click the button right next to the drop down menu to show the Network Analyst Window. You’ll be using it a lot.
The window will have a list of Stops, Routes and Barriers, all empty for now. The first step is to populate your list of stops. You can create stops manually by activating the Create Network Location Tool and clicking on the screen while your pointer is in the flag with crosshairs shape. Or, you can load them from another layer. That’s what we’re going to do. Select the 12 sites that you see on screen. (This step is important because you don’t want to be working with all the recreational areas in the entire USA. I skipped it the first time and proceeded to crash my computer). Once you’ve selected your subset, right click on Stops and select Load Locations…
It’ll grab their name from a name field if one exists, and nicely number them for you. Now click the Solve button. It’s generated a topsy-turvey looking route, due in part I’m sure to all the one-way streets downtown.
If you click the driving directions button, you’ll see the route has 80 steps and takes about 14 minutes.
We can probably do better than this. There’s an obscure button that will allow you to change some of the parameters used to calculate your route. It’s hidden on the top right of the Network Analyst Window.
The Analysis Settings tab has a Reorder Stops to Find Optimal Route box that is checked off by default. This means that if you added the stops manually, it will visit them in the order you created them, which might make sense. But if you loaded them from a file like we did, it will visit them in the order of random numbers it assigned, which definitely doesn’t make sense. Check the Reorder Stops box on, and check off the two boxes under it: Preserve First and Last Stop.
Now, the software will have the freedom to plan an itinerary that maximizes efficiency. Let’s see what difference it makes. After you’ve pressed OK to dismiss the Route Properties box, press the Solve button again.
That’s much better! It now wants me to visit the Children’s Museum first instead of last. And it has rearranged the entire middle so it’s not a mess. This new route saves us 5 steps, 3 miles and 4 minutes.
Service Area Tool
Now let’s pretend we work for Woodland Park Hospital and we want to know the area around it than an ambulance can reach within 5 minutes. We would select New Service Area from the drop down menu, and then load that one hospital point into the Facilities list.
I personally like to go into the Service Area Properties box and change the Polygon Type to Detailed (under the Polygon Analysis tab) because, why wouldn’t you? It’s set to calculate for 5 minutes by default, but you can change this to anything by editing the value in the Default Breaks box (under the Analysis Settings tab). If you want to solve for multiple times, separate them by commas. Here’s how it looks if I solve for 1, 3, and 5 minutes with high detail:
By the way, if you want to change these colors you can edit the symbology just like you would any other layer in the table of contents. And, if you wanted to, you could turn off the OneWay restriction because ambulances don’t need to obey it.
Closest Facility Tool
For the final scenario we’ll pretend to be home owners in the southwest portion of the city. We are planning a family and would like to determine the closest school to our house. The first step to solve this problem is to select New Closest Facility from the Network Analyst drop down menu, then press the Create Network Location tool button, and click on the map at the location of the house.
The new location will automatically be added to the Facilities list, but you want it in the Incidents list. Either drag it down there, or make sure the Incidents list is highlighted before you use the Create Network Location tool. The point will be automatically named Graphic Pick 1. You can double click on its entry in the Incidents list to rename it.
I purposefully picked a location that seems equally close to three different schools. The tool will tell me which is actually the fastest to drive to, not the closest as the crow flies. I selected a bunch of schools around it (more than just those three, to be safe) and loaded them into the Facilities list. Then I pressed the Solve button. The result is that Markam School is the closest at 0.7 miles away.
The tool shows you the route to get there, too. You can open up the Directions Window to see driving directions as well.
These are pretty handy, aren’t they? I’m in the mood to plan a vacation now just so I can play with the Route Tool some more!
Tags: driving directions, network analyst, service area
Creating a network dataset
This is post #2 in a three part series. It is meant to be read after Preparing your Data to Become a Network Dataset.
A network dataset can be created from a shapefile or a feature class. If you are creating one from a feature class, it needs to be in its own feature dataset. Make sure you have the Network Analyst extension checked on, then right click on your streetline shapefile or the feature dataset that holds your streetline feature class, and choose “New… Network Dataset”.
After you give it a name, you will have the option the set Connectivity. Select the connectivity policy that you decided on before.
If you have designated z-levels the way I described, the next screen you encounter will automatically have these elevation settings. Keep them.
I haven’t tackled turn modeling yet, so select “No” on the next screen unless you’ve figured that out. (In which case, do you want to do a guest post?) 
Now you are at the Attributes screen. If you’ve designated one-way streets the way I described, they will already be set up as a Oneway Restriction attribute. In addition, you’ll have a RoadClass Descriptor Attribute. Network Analyst uses this attribute to make the wording for turn-by-directions read correctly. There are 5 categories:
1 Local Street
2 Highway
3 Ramp
4 Ferry
5 Roundabout
It is automatically configured to try to inspect your attribute table and pull out these categories, but there’s no guarantee it’s going to work right. You should look at the code to check. Select it and click on the Evaluators button. Then select the first one in the list and click on the Evaluator Properties button.
Here’s what appears in my Pre-Logic VB Script Code box:
rc = 1 'Local road
Select Case UCase([CFCC])
Case "A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18", "A19"
rc = 2 'Highway
Case "A60", "A63"
rc = 3 'Ramp
Case "A65", "A66", "A68", "A69"
rc = 4 'Ferry
Case "A62"
rc = 5 'Roundabout
End Select
It is picking attributes out of my Census Feature Class Code (CFCC) column. But is it not classifying them in the same way my metadata says they should be. I have some evidence my metadata is out of date, so maybe the code knows better, but I can’t be sure. So, I replaced the code with this, which I am sure about:
If [FTYPE] = "FWY" Or [FTYPE] = "HWY" Then
rc = 2 'Highway
ElseIf [FTYPE] = "RAMP" Then
rc = 3 'Ramp
Else
rc = 1 'Local road
End If
I didn’t have anything that was obviously ferries or roundabouts. Make sure rc, or whatever variable you use, is in the “Value =” box below the code window. Make the same changes for the To-From direction as well.
Now it is time to create your Cost Attributes. Press the Add… button and create a Distance Attribute. Make sure the Usage Type is Cost and that you set the units and data type to match your length field. Uncheck the Use by Default box unless you don’t have information to create a Time attribute.
It will have a yellow exclamation point beside it until you press the Evaluators button and define how it should be calculated. Set the Type to Field and the Value to your length field for both directions.
Follow the same procedure to create a Time attribute, except this time set the units to Minutes and keep Use by Default checked on. When you get to the Evaluators screen, set the Type to Field and then click the Evaluator Properties button so you can enter an expression (unless there is already a time field in your attribute table). I used the expression ( [LENGTH] / [Speed] ) * 60 in the “Value =” box. This works if you have a field for length in miles, speed in miles per hour, and want time in minutes. Adjust accordingly.
Next create a Hierarchy attribute (Usage Type: Hierarchy). You will probably want to leave the Use by Default box checked, but I will leave that up to you. (These default settings can still be changed when you run scenarios later). If you created a new field for it, then calculate it using that. Otherwise you’ll have to code a classification expression like I did for the RoadClass attribute… only this time with categories 1-3.
Finally, if you thought of any potential restrictions to travel earlier, create Restriction attributes for them now. Create a different attribute for each type of Restriction. I set one: “Private”, for private roads and driveways, so I could avoid them if need be (like if I was driving a commercial vehicle). I didn’t use it by default because most of the time this won’t be an issue. Restriction attributes are based on a Boolean value (True or False) so my classification expression looks a little different:
Here’s my list of attributes in the end:
If you’ve done everything right thus far, you’ll have the option to configure driving directions. (It didn’t let me the first time around, I think because I didn’t create a Distance attribute). These are cool and you definitely want them! Click the Directions… button to see/change the settings. It’s pretty smart and gets everything right for me automatically. You can click anywhere to bring up drop-down menus and change things.
I didn’t have anything in my attribute table that helps with labeling highway shields (like route numbers). But I do have a field that says which city/county each road segment is in. So I pointed to that as the boundary field on the Boundary tab. Now, the driving directions will tell me if I cross over from one jurisdiction to another.
You’ll get a summary at the end, and then it’s time to click FINISH! Go ahead and Build when it prompts you. Then, add your new network dataset to ArcMap because it is time to play.
Stay tuned for installment #3: Doing network analysis
Tags: driving directions, junction, network analyst, VBA
Preparing street data to become a network dataset
Problem: You are unfamiliar with the Network Analyst extension and you would like to maximize the potential of your streetline data.
Recently I encountered my first network dataset. This is a special file created by, and meant to work with, the Network Analyst extension.
Route planning is the first thing that pops into my mind when I think about what a GIS does, yet I’d never even turned on that extension before. This is going to be a very basic introduction for people who might be in the same boat: interested in what Network Analyst can do, but don’t know where to start.
I am going to post it in three installments
1. Preparing your data to become a network dataset
2. Creating the network dataset
3. Doing network analysis
The first thing you will need is a streetline layer. The more information it has in it, the better. It should at least have road names, road types, and direction of travel (for one-way streets) designated, or its use will be very limited. You can use the ESRI StreetMap data if you have access to it, which is optimized for use with Network Analyst and even comes with a pre-made network dataset. If you’re going to do that you can skip ahead to installment #3. Otherwise, read on:
Preparing your data to become a network dataset
I downloaded a streetline shapefile from the City of Portland GIS Data Catalog.
The fundamental bit of information network analyst needs in order to work is a Cost Attribute. Cost is either the travel time or distance on each road segment. It is the thing the model is trying to minimize when calculating the best route. Distance is easy to add to any dataset. Create a field (type: double) if one is not there already, and then use the Calculate Geometry tool. Chose miles or kilometers, whichever you prefer, as your length unit. Note: your dataset needs to be in a projected coordinate system for this to work.
If you know, or have a basis upon which to estimate, the speed of travel on each road segment, you can also create a Cost Attribute based on time. This will make your results much more meaningful. My dataset designated road types, which allowed me to make an educated guess about speed limits. I added a Speed field and populated it with my miles per hour estimate. The Insurance Institute for Highway Safety has a table of maximum speed limits in each state. Oregon’s is 55 for urban interstates. So, I used 55 for highways and freeways, 45 for major/primary roads, 35 for secondary roads, 25 for minor/local roads and ramps, and 15 for driveways.
These speed and distance fields in the attribute table will be used to calculate a time attribute later in the process. Moving on:
The next fundamental attribute to have is Hierarchy. Network analyst uses it to create more realistic routing based on the fact that people prefer to travel on larger roads. It is an integer value that ranges from 1 (major) to 3 (minor) based on road class / road level. If your roads are already labeled as primary, secondary, and tertiary, you’re set; otherwise, you’ll have to make some classification decisions. You can add a new integer field to the attribute table for it (that’s what I did). Or, you can calculate it later, like I’ll explain how to do with the time attribute.
Next, you will need to flag one-way roads with their direction of travel. I struggled with this for quite awhile before stumbling upon a procedure that works correctly, so I recommend doing it just this way:
- Create a new string field in the attribute Table called OneWay.
- Change the symbol of your roads to a cartographic line with the arrow pointing right. This will allow you to see what direction the line was digitized in.
- Bring up a Google map of your area, or something else that allows you to see the direction of travel.
- Put an “F” in the OneWay column for any road that has been digitized along with the direction of travel.
- Put a “T” in the OneWay column for any road that has been digitized against the direction of travel.
- Use an “N” for any roads that don’t permit travel in any direction.
- Use any other value for 2-way streets. (I used “X”).
(click on image to see full size)
When you create your network dataset later, Network Analyst will recognize this formatting and will know what to do with it.
Finally, examine the way your road lines were digitized, especially around overpasses and tunnels. If they have a unique segment between each intersection (that is, they break anywhere a person could theoretically turn), then you can use End Point connectivity—the default. If, however, your roads are long lines that don’t break when they cross, you will have to use Any Vertex connectivity. The type of connectivity you chose tells Network Analyst where to create Junctions. Junctions are a special point layer that will be generated along with your network dataset. They are created wherever a person could transfer from road to another.
Regardless of what connectivity you chose, you will need a way to designate overpasses and tunnels so that people will not be routed to turn onto roads that are inaccessible because they pass below or above them. This is done using z-levels. Create two new integer fields to hold your z-level information. Call one “F_ZLEV” and the other “T_ZLEV”. The field prefixes stand for “From” and “To,” respectively. First, designate all ground level roads as z-level 0 in both fields. Then designate overpasses as 1 and tunnels as -1 in both fields. (Hopefully, your attribute table will tell you where these are. Mine had a “Structure Type” field). Now, symbolize with a cartographic line again as I described in Step #2, so you can see the digitization direction. Change F_ZLEV to 0 at the beginning of raised or lowered segments where they meet level road. Change T_ZLEV to 0 at the end of raised or lowered segments where they meet level road. The idea is that a route will only go through a junction where the z-levels match.
(click on image to see full size)
You are almost ready to move on to the next step, but before you do, think about any other potential restrictions to travel. Do you have roads that are marked as “under construction,” toll roads or private roads? Make a note of how they are designated in the attribute table so you can flag them later.
Stay tuned for installment #2: Creating a network dataset
Tags: cost attribute, network analyst, one way street, z value