MapWindow Measuring Tool and curious NAD 83 State Plane Coordinate System inconsistencies in .prj files

I have been using MapWindow recently, and it seems like a very good program, especially considering it’s free. As my employer no longer provides a viewer to our customers, we have been recommending MapWindow…but…there is a problem. Many of our customers use the California State Plane Coordinate System with units of U.S. feet and, for whatever reason, the measuring tool in the current version of MapWindow (4.8.6) doesn’t seem to be able to deal with feet. I started noodling around with shapefile projections in ArcMap, Global Mapper, and MapWindow and checking out the .prj files in Notepad to see if I could find any clues.

So, I created a shapefile and made a copy. I defined the projection of one in ArcMap, and of the other in MapWindow. Arc was able to display them both correctly, but did not recognize that the projection is actually the same. When I pulled both into MapWindow, the ESRI one failed to reproject and was not displayed. When I opened both in Global Mapper, no problem. So then I exported a new one from Global Mapper, and tried opening all three simultaneously in MapWindow. MapWindow displayed the Global Mapper and MapWindow ones, but again the ESRI one failed. All three open fine simultaneously in both ArcMap and Global Mapper, but the addition of the MapWindow defined layer triggers a warning that it has a different geographic coordinate system than the other two.

A peek at the projection information illuminates that there is some inconsistency in how even the datum is described:

From MapWindow:

NAD83 / California zone 5 (ftUS)
Projection: Lambert_Conformal_Conic
false_easting: 6561666.667000
false_northing: 1640416.667000
central_meridian: -118.000000
standard_parallel_1: 35.466667
standard_parallel_2: 34.033333
latitude_of_origin: 33.500000

Linear Unit: US survey foot
NAD83Datum: North_American_Datum_1983

From ArcGIS:

Projection: Lambert_Conformal_Conic
False_Easting: 6561666.666667
False_Northing: 1640416.666667
Central_Meridian: -118.000000
Standard_Parallel_1: 34.033333
Standard_Parallel_2: 35.466667
Latitude_Of_Origin: 33.500000

Linear Unit: Foot_US
GCS_North_American_1983Datum: D_North_American_1983

From Global Mapper:

Projection: Lambert_Conformal_Conic
false_easting: 6561666.666667
false_northing: 1640416.666667
central_meridian: -118.000000
standard_parallel_1: 34.033333
standard_parallel_2: 35.466667
scale_factor: 1.000000
latitude_of_origin: 33.500000
Linear Unit: Foot_US

GCS_Geographic Coordinate System

Shouldn’t this kind of thing be standardized?

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)

Cell counts and flow direction, part 2

Let’s see…where did I leave off?

After posting a note on facebook hypothesizing that perhaps the graph shown in the previous post showed that the flow direction algorithm was biased toward the ordinal directions.  Brian agreed that there was an artifact here, perhaps as a result of the algorithm, or possibly of the gridded nature of the elevation raster itself.  Brian suggested a couple of tests, involving conversion from raster to point and reinterpolating the raster, converting to contours and then using topo to raster, or smoothing using focal statistics.  I did not try any of these, because I did not want to open the interpolation can of worms here, per se, because I thought it might essentially introduce another variable and muddy the source(s) of bias further.  I had reprojected the data from NAD83 to State Plane, because otherwise slope never comes out right, and this projection introduced an artifact mesh.  After a little back and forth on this, Greg chimed in that DEMs are highly inaccurate to begin with as a result of data collection.  , and that conversion back and forth only propagates the error.

I realized that the sine wave implies a D-infinity flow algorithms–that is, flow could be calculated in any direction from the cell center.  This is, of course, not what ArcGIS actually does.  What the graph should really show is 8 columns to represent the 8 discrete directional possibilities from one cell to any of its eight neighbors.  What occurred to me at this point is that for any of these 8 possible directions, the cell count will not be the same at a given distance.  If, for example, one measures the distance from the center of a cell to the center of each neighboring cell in the ordinal directions (diagonally), then the distance will be the cell width multiplied by the square root of two, and we can call the cell count 1.  Assuming that cells can be partially counted, then the same distance in the cardinal directions will cover about 1.414 cells.  Therefore, the pattern I noted may be nothing more than the geometry of a square, causing cell counts to be relatively lower in the ordinal directions for exactly the same flow length.  I thought a way to test this might be to clip a circular extent and rotate it.  It should be possible to remove the effect of this geometry by using a factor of the square root of two, and then to see whether the same number of cells is calculated for each angle of rotation (corresponding to the 8 directions in question) when the angle is accounted for.  I made a graph that shows how I believe the pattern should change for an identical area as it is rotated through a full rotation at 15-degree increments.


Also, I made a couple short animations in Windows Live Movie Maker showing the flow directions as the circular clip from the DEM is rotated.  These aren’t so much for analytical purposes as for fun.  They were more fun with music, but I removed it to make the files smaller:




Cell counts and flow direction

Recently, as I was playing around with the DEM for what I will undoubtedly inadvertently refer to as “my watershed”, as if I own it. The watershed in question is a small one…I don’t know the area off the top of my head but on my other screen here I have ArcMap open, so…looks like 707392.5 square meters, according to the Shape_Area field of the watershed polygon I converted from the basin computed from the flow direction based on the filled 1/3″ NED. Bit of false precision there, I’d guess…for convenience and because really what I wanted to express is the order of magnitude, let’s just call it 707000 sq m, which I guess is what… 0.7 square km, as a ballpark figure. It’s at the upper part of the Laguna Canyon watershed, in the sort of northeastern San Joaquin Hills, conceptually defined as the area draining to Barbara’s Lake.
I was looking at the DEM for a squarish area covering an area several times larger than the watershed. Again, I have Arc open, so…about 800 cells on a side, so about 8 km on an edge, 64 sq km overall. And I was looking at the flow direction raster for the area of the watershed and thinking that because the bulk of drainage, that is, the accumulated overland flow, should go southwest toward Laguna Canyon and out to the Pacific, to the adjacent parallel canyon and in the same direction, or to the northeast corner of the map where it goes out onto the plain and turns south, ultimately to the Pacific as well. Clipped in tight, it should all go southwest. But how does the flow accumulation, I wondered, compare with the flow direction. So I opened up the attribute tables and took a look at the count of each cell classification in the flow direction raster. I won’t go deep into the various tools I mention here…rather, I would refer you to ESRI’s online desktop help page, where you can find more detail. I copied the numbers over to Excel and made first something that looked kind of like a scree plot, for those who know what that is. The labels below the horizontal axis were ordered highest count to lowest. To my initial surprise, the majority of cells were classified east, followed by west. The other directions all had much lower values, but north and south were next. I could see most of the flow being east and west, perpendicular to the ridges and the channels–that makes perfect sense, but why then north and south, rather than northeast, southeast, northwest, southwest? I sorted in numerical order so that the axis began with due east and went around in a circle clockwise. With a smooth line fit through the numbers, it looks a bit like a sin wave with some kind of amplitude boost every other upward peak. Really, because ArcGIS calculates flow to the 8 neighboring cells, a more appropriate graph would reflect the discrete nature of the possible values, but I include the graph here uncorrected in this regard.
I posted a note on facebook showing this graph, and pondering whether this pattern, each peak corresponding to one of the cardinal directions (E, S, W, N) and each trough one of the ordinal (SE, SW, NW, NE), was an artifact of the gridded data, perhaps bias introduced by something in the flow direction algorithm. I received some feedback from colleagues at Long Beach, Brian Nagy and Greg Ziolkowski, stimulating and combining with some of my own evolving thoughts on the matter, which I well get into next blog.

An introduction

The little “About Me” gadget over there has a little blurb about me…grad student, geographer, blah blah blah. To the academics and professionals among you who find my casual writing style shocking (am I daring to address you, as though I am the Abby behind Dear Abby?), let me plainly say that I will be keeping this rather informal. The author wishes to avoid the detached narrative voice of the third person (omniscient?) who so often describes what methods were carried out (without an actor? of their own agency?) or what results imply (free from the presence of an interpreter?). Forgive me. May I point out that in a sense, I am first a writer, then a geographer, as my first degree is a B.A. in English, writing minor. I was trained to despise passive voice, but it so often seems the only way to avoid “I” or “we” or “me”. It was my years of wondering what to write about that drove me back to college–partially, anyway–and only then did I stumble into the discipline of geography. I will write this blog informally, and stream of consciousness (punctuation be damned!), as it is my blog, not an article for submission to a peer-reviewed journal.
I would like to expand upon my rather expansive list of interests listed over there, “About Me”. I do consider myself a physical and environmental geographer, although I would like to stress that this does not mean I have zero interest in the human, the cultural, the social. I just don’t know as much about it, because there isn’t enough time in the day to learn all that I want to learn, and so I must prioritize. I prioritized geography over geology, though I am always eager to learn more about rocks and minerals and geochemical cycles and the force of flowing water in shaping the landscape, both in terms of erosion and deposition, dissolution, transport, and sedimentation and in terms of the relationship with life, with plants and animals and humans. I have a fondness for bugs…insects, arachnids, and the bugs of the sea…invertebrates in the tidal zones…but unfortunately I have not a biology class since high school. I am certainly capable of learning things on my own, but again, not enough time in the day to learn everything. Bugs and little invertebrates, microbes, and so on are so nearly invisible, or thought creepy, but they are little ecological engines playing such crucial roles. I have two cockroaches…Gromphadorhina portentosa. But back to geography…
So similarly, I have learned a bit here and there about water policy, watershed management, land use, and environmental justice, political ecology, and the like. I have learned a bit about physical hydrology and hydrogeology. And yes, I have spent quite a bit of time with GIS, at this point. Initially, I took a GIS class at Santa Monica College simply to make myself more marketable for environmental jobs. My sister, an anthropologist, had studied environmental science and environmental studies, and reported to me that she and her fellow students in those majors were being asked about GIS skills when seeking jobs or internships. So I figured I had better get on board. In that first GIS class I immediately understood…a geographic information system is a powerful tool for environmental analysis. It can be used as an analytical tool in itself, or be used for pre-processing data for other models, and for assembling output maps for presentation.
Fast forward, and I have just finished the last class I needed for a GIS certificate. I have also finished all of the coursework required for the M.A. (and then some)–all I have left is the thesis. My thesis topic is DEMs as a source of uncertainty in modeling the Barbara’s Lake watershed. This is an area in which my advisor and a colleague in Geology have conducted research, and in which they hope to continue to conduct research. I hope to find out how the uncertainty in the digital elevation model impacts the actual predictions of runoff patterns in the watershed and pollutant loadings in the lake.
I suspect it will be negligible compared to other sources, considered strictly in terms of accuracy, that is, whether the elevation one would actually find at a given location in the field is, within the error of the measuring instrument and method, what is shown on the map. The most obvious error encountered so far actually appeared in work a fellow student, Brian Nagy, has been doing in creating a digital model to serve as the guide for a solid terrain model. The DEM, it seems, was created from a USGS topographic map dating to the 1970s. In the time since the source map was created, cuts and fills and leveling have changed the lanscape. In terms of real world application at this site, the temporal sampling interval is the limiting factor on accuracy, here, not the spatial resolution or the intrument precision or what have you. If the landscape was unchanged from the source map, some uncertainty would still exist from these other sources…this is the uncertainty I wish to measure the impact of. I’ll let you know how it turns out.