One of the most exciting new features in ArcGIS 10 is the introduction of iterating tools in ModelBuilder. The ability to step through a list has always been a mainstay of programming and computer science, and has been readily accessible to any GIS programmer with a little Python or C# experience. In previous versions of ModelBuilder, iteration was possible through the clever use of count variables and batch processing. But now it is available graphically on the desktop, and can be added to a model with a single click.
To illustrate, consider the following example. We begin with 2 files, a set of building footprints (derived from LIDAR), and a collection of census block groups. In order to determine the built-up density of each block, we need to determine which of the building footprints are in each block, and then calculate the total building area. To accomplish this, we need to step through each block, extract (Select) it from the dataset, and then use it to Clip the building footprints. Then we can use Summary Statistics to sum the building area. I’d also like to add the block group ID as a field in the statistics file, for later reference, and this can be done with calls to Add Field and Calculate Field.
Although it sounds difficult on paper, this is a fairly typical task for ModelBuilder. There are only 5 actions, and the workflow is largely linear. The trouble is that we want to do this 516 times! This is where iterations become handy. One of the iteration methods, “Iterate Field Values” allows you to step over each row in a field and pass the value of the field as a variable. In my example, the variable holding the block group ID is simply called Value, and I can reference it anywhere in my model by calling %Value%:
Click the image to view a larger version
Notice in the above graphic that the iteration is not attached to the main workflow. Instead, it appears to run parallel. Although this notation may seem bizarre to seasoned programmers, I assure you that it is identical to the behavior of a typical for-each loop.
At this time, only one iteration can be used for each model. For more information about how to use iterations, check out the online help.
TriMet, the Tri-County Metropolitan Transportation District, serving Washington, Multnomah, and Clackamas counties in Oregon, is a leader in open source technology. They still use ESRI desktop applications, but they are using open technologies, including PostGIS and OpenLayers, to power their web services.
Article:
http://www.govtech.com/e-government/Open-Source-Software-Oregon-Transportation.html
If you’ve ever needed to export a shapefile or feature class to a GeoJSON text file, there’s a handy ESRI-supported Python script available. It’s easy to install into your ArcToolbox, and works quite well:
http://arcscripts.esri.com/details.asp?dbid=15545
So what is GeoJSON? JSON stands for JavaScript Object Notation, an encoding format for web object data. GeoJSON is just a JSON object that contains spatial information. In a sense, it’s a highly portable web container for geographic data. Many web mapping programs can parse it, including PostGIS and OpenLayers. You can also use it when programming with the Google Maps API.
Dr. Matthew Wilson, a Geography professor at Ball State University (Indiana) recently posted a reading list for critical GIS. What exactly is critical GIS? It’s not a synonym for emergency GIS. What Dr. Wilson actually means is critical thinking about GIS – not how to use GIS, but why to use GIS. The list is broken up into a timeline, so you can trace the evolution of critical GIS through the readings. I admit that I’ve only read a handful of the selections, but I can attest to the value of this line of study. Some of the pieces are in journals that you or your organization may have on hand already, and most of them can be found online or in an academic library. It’s never too early to work on that summer reading list!
http://criticalgis.blogspot.com/p/critical-gis-bibliography.html
The hypsometric contour interval on the classic USGS quadrangle maps is 40 feet. Although it’s not a nice round number, I think it looks fantastic at that scale. It’s good at representing many different terrain types. Which is why USGS contours have become a bit of a standard for many GIS and cartographic applications.
But what if you are mapping at a much larger scale? For a county-sized area, 200 feet contours may be more appropriate. Unfortunately, it’s not always trivial to recontour the original surface, especially if the data is not readily available.
The good news is that you can easily extract the 200 feet contours from the 40 feet contour dataset. The key to this process is Select by Attributes in ArcGIS. A contour dataset is really just a set of polygons with an elevation field (named ‘ELEVFT’ in my sample dataset). Selecting only the contours that are evenly divisible by 200 requires the modulo function. The modulo function takes two numbers, divides them, and then returns the remainder. If the remainder is zero, that means the two numbers divide evenly. In ArcGIS SQL syntax, this is what your selection statement looks like:
MOD(“ELEVFT”,200) = 0
Once the contours have been selected, it’s easy to create a new shapefile and add it to the map. You can also extract the 1000 feet contours, and plot them as index contours. I personally like to use 1pt black lines for contours, and 1.5pt black lines for index contours. On a busy map, I also like to reduce the contour layers to 60% transparency.
One source of data that I feel is under-utilized is the USGS Board of Geographic Names. The BGN is a (mostly) comprehensive list of definitive feature names, covering most of the United States. The features are organized by state, by county, and by feature type, and each one has a latitude and a longitude.
The only drawback to the BGN is that the files are pipe-delimited (the ‘|’ character is called ‘pipe’) text files. For example, here’s the first five entries from the Hawaii file:
247074|Pacific Ocean|Sea|CA|06|Mendocino|045|391837N|1235041W|
39.3102778|-123.8447222|||||0|0|Mendocino|01/19/1981|06/21/2010
358293|Pearl and Hermes Atoll|Island|HI|15|Honolulu|003|275000N|1755000W|
27.8333333|-175.8333333|||||0|0|Unknown|09/30/2003|
358294|Laysan Island|Island|HI|15|Honolulu|003|254615N|1714415W|
25.7708333|-171.7375|||||3|10|Unknown|09/30/2003|
358295|Barking Sands|Beach|HI|15|Kauai|007|220418N|1594652W|
22.0716667|-159.7811111|||||0|0|Kekaha|02/06/1981|
(Yes, ladies and gentlemen, the Pacific Ocean is in Mendocino County, California.)
This format, while containing a lot of great information, is not exactly ArcGIS friendly. Fortunately, it is really easy to transform this file into a comma-separated file, which can easily be manipulated in a spreadsheet, using awk. Since there are 20 columns in each entry (whether the column is populated or not), the awk script will look something like this:
awk ‘BEGIN {FS = “|”; OFS = “,”} {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16,$17,$18,$19,$20}’ HI_Features_20101203.txt > HIfeatures-awk.txt
That’s all there is to it. The file can now be loaded into a database and plotted on a map.
If you ever get the opportunity to work with bathymetric data, you may encounter CARIS files. CARIS is an enterprise GIS suite that is well suited to marine and hydrographic applications. Many GIS offices don’t have access to CARIS, but that doesn’t mean the data is out of your reach. CARIS offers a free data viewer which you can use to export the data to a more versatile format – ASCII text.
The first thing you’ll want to do is download the CARIS EasyView program. You’ll need to provide name and contact information, and you cannot download the program until you’ve received a confirmation code via email. Installation is typical for Windows programs, and opening the CARIS file (which will probably have a .car or .car0 extension) is fairly straight forward. The first thing to do in EasyView is select File > Export Surface to ASCII, which brings up the following dialogue box:
Figure 1: CARIS EasyView export dialogue
If there are multiple surfaces in your file, you’ll have to export them one at a time, by selecting them in the dialogue box. In my file, there was only one surface to convert, so picking the right one was easy. You’ll want to set the position units to Geographic DD – decimal degrees. Although my data had a number of different features which I could output, I was only really interested in the depth. Include any features in the output file by adding them from the list on the left to the list on the right. One of the key attributes to choose here is the delimiter. I chose the comma, making the output file a de facto CSV – which is readable by both ArcGIS and major spreadsheet programs. My research indicates that the final attribute, Z-axis convention, has no effect on the data output – but it’s probably best to stick with “Up is positive” just in case.
The textfile created by this process is a comma delimited set of coordinates.
As an example, from the file I was working with:
44.031723N,076.495925W,8.872
44.031723N,076.495925W,8.817
Unfortunately, these numbers are not very useful to us, as ArcGIS uses negative numbers to represent longitude west. The combination of number and letter in the latitude and longitude fields will make projecting and manipulation difficult. And it’s better to have negative depth values (and positive altitude values). So I’m going to use awk, a command-line tool, to process the file:
awk ‘BEGIN {FS = “,”; OFS = “,”} {print $1*1,$2*-1,$3*-1}’ carisexport.txt > carisawk.txt
Without getting into too much detail, this multiplies the first column by 1, and the second and third columns by -1, using a comma as the field separator. (Multiplying an alphanumeric string by a number in awk recasts the variable as a number.) Now we need to add simple headers to the file, signifying names for each field. Open the file in a text editor – Windows Notepad will probably reject the file for being too long, so a more functional text editor (such as EditPad Lite) is required. Latitude, longitude, and depth are good names for the fields. Here’s the first 3 lines after running awk and adding headers:
latitude,longitude,depth
44.031723,-76.495925,-8.872
44.031723,-76.495925,-8.817
Now the file can be opened using Tools > Add XY Data in ArcGIS. Select latitude as the Y Field and longitude as the X Field, and don’t forget to define a projection. Export the layer as a shapefile, and it can be plotted, rasterized, or spatially analyzed as you see fit.
If anyone from CARIS is reading this, it would be very helpful if you added simple headers to the ASCII export function, and used standardized representations for latitude and longitude values. Maximum compatibility with companion tools is a strong selling point for any proprietary software.
As a spatial professional, projections and geodesy have become second nature to me – often to a point that I take them for granted. In many cases, the projection becomes like a font in a presentation – if you used Helvetica on page one, you had better use it on page two. Generally, the consistency is more important than the overall effect. When it comes to fonts, there’s only a handful of professional looking typefaces to choose from, and one is little different from another. Projections can sometimes feel this way, especially if you spend most of your time looking at similar projections, such as adjacent state plane systems.
Occasionally, though, you get a chance to see just how much difference a projection can make. Over the weekend, I downloaded a grid of bathymetry and hypsometry of Michigan’s Great Lakes from NOAA’s NGDC Grid Translator. The grid was in geographic form, but had neither projection nor datum applied to it – just raw latitude and longitude. The raster I created from it looked somehow flattened and lazy:
Figure 1: The Great Lakes, unprojected
This is a reasonably accurate view of the area, but it’s not the view we’re accustomed to seeing. Although a worldwide coordinate system might make sense when viewing a whole hemisphere, a localized projection that minimizes certain kinds of distortions is more suited to a smaller area. As it happens, most of Michigan’s Great Lakes happen to fall within Universal Transverse Mercator (UTM) Zone 16. Which gives a much nicer looking image:
Figure 2: The Great Lakes, projected in UTM16
This is more like the Great Lakes most of us recognize. The UTM projection preserves the shapes and angles our eyes have grown accustomed to, at the expense of distance and area. Notice, of course, how this image looks more like a slice from the outside of a sphere, unlike a typical rectangular map. Although I wouldn’t normally use this kind of shape for a presentation map, I think it illustrates nicely just how the UTM projection modifies the mathematics of the raster to create the correct image.
PostGIS is a wonderful suite of tools, which really anchors the open source web mapping toolkit. However, it doesn’t always integrate well with an ArcGIS desktop environment. One obvious source of conflict relates to projections – most of ArcGIS’s projections are infinitesimally different than their real world counterparts. As an example, consider the Washington North State Plane coordinate system. Although the official unit of the State Plane system is the meter, the projection has also been defined using the foot. The Spatial Reference ID (SRID) of the meter version is 2855, while the foot version is 2926. However, the version of the coordinate system that ArcGIS uses is SRID 102748. Unfortunately, this projection, which almost all Washington government agencies use, is undefined in PostGIS.
The good news is that you can add new projections to PostGIS. It’s neither simple nor intuitive, but it is quite possible. The following is a set of instructions for adding projection 102748. I’m assuming you have a default installation of PostGIS and total access to your server.
Adding a projection is accomplished by modifying the spatial_ref_sys table, which is part of the template_postgis database. But by default, that database doesn’t allow connections. So you’ll need to connect to template1, as the user postgres:
psql -d template1
Modify template_postgis so that it allows connections:
UPDATE pg_database SET datallowconn = TRUE WHERE datname = ‘template_postgis’;
Now you can connect (switch) to template_postgis:
\c template_postgis
The code to actually add the projection is long and complicated. Fortunately, spatialreference.org supplies us with the perfect code. Just click ‘PostGIS spatial_ref_sys INSERT statement’ and copy the code:
INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 9102748, ‘esri’, 102748, ‘+proj=lcc +lat_1=47.5 +lat_2=48.73333333333333 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs ‘, ‘PROJCS[“NAD_1983_StatePlane_Washington_North_FIPS_4601_Feet”,
GEOGCS[“GCS_North_American_1983”,DATUM[“North_American_Datum_1983”,
SPHEROID[“GRS_1980”,6378137,298.257222101]],PRIMEM[“Greenwich”,0],
UNIT[“Degree”,0.017453292519943295]],PROJECTION[“Lambert_Conformal_Conic_2SP”],
PARAMETER[“False_Easting”,1640416.666666667],PARAMETER[“False_Northing”,0],
PARAMETER[“Central_Meridian”,-120.8333333333333],PARAMETER[“Standard_Parallel_1”,47.5],
PARAMETER[“Standard_Parallel_2”,48.73333333333333],PARAMETER[“Latitude_Of_Origin”,47],
UNIT[“Foot_US”,0.30480060960121924],AUTHORITY[“EPSG”,”102748″]]’);
All done, let’s switch back to template1 and reset template_postgis:
\c template1
UPDATE pg_database SET datallowconn = FALSE WHERE datname = ‘template_postgis’;
These changes should take effect once the server is restarted.
So is it worth all the work? Sorta. Using one of the more standard SRIDs is plenty accurate – I think the only real difference is a few decimal places. But I think it’s important to be vigilant, especially when dealing with projections. Precision, accuracy, and internal consistency are the cornerstones of science and repeatability.
In my previous post, “When Scale Bars Lie”, I mentioned the dangerous things that can happen when the map unit is not set correctly. As I explained, it can lead to abnormal values, which are often manifested in scale bars. This week, a friend of mine was using GIS to calculate the areas of some urban parcels, and getting values in the billions. Once again, the map units were to blame. Setting the map units and display units solved his problem instantly.
As I mentioned before, the easiest way to prevent this issue is to be mindful of projections. I’ve been burned by projection problems before, and on one occasion I had no choice but to restart the project from scratch. So now, I am extraordinarily mindful of projections. Before I begin a project, I choose which projection to use (or ask the clients which projection they prefer). Once a projection has been decided upon, it becomes law – every piece of data is reprojected before being inserted into the project. Admittedly, this is overkill, because ArcGIS will throw up an alert and offer to fix major projection conflicts, and minor conflicts usually don’t add up to much. But I’ve decided that it’s better to be safe than sorry – especially when dealing with unprojected datasets. As the old adage goes, measure twice, and cut once.
