Geodesic Distance in ArcMap

One of the things I like about ESRI ArcMap is the ease with which you can measure features in your data. Recently I wanted a a simple way to measure geodesic distances of some projected data. The measure tool in ArcMap, by default, measures distances in projected units in projected data, and geodesic distances only for data displayed in geographic lat/long. If you hold down the shift key while measuring a distance, Arc will calculate geodesic distances regardless of projection. What I really wanted to be able to do, however, was to be able to draw lines in a shapefile, and then have Arc calculate the geodesic distance of those lines. This is something that Arc cannot do, without you first unprojecting (or more correctly, reprojecting) your data to geographic coordinates.

Luckily, I found this: http://137.227.239.67/pigwad/tutorials/scripts/

Looks like some of the folks at ESRI got together with the USGS a designed a plug-in that will let you, among other things, calculate geodesic distances from projected data. It uses the parameters specified by your projection to unproject the data into geographic coordinates and then calculates the geodesic distance on the fly. Since it uses the semi-major and semi-minor radii that you specify in your projection, the distortion should be minimized. This tool is definitely a time saver, as it does the reprojection for you only on the data you are calculating, and enables you to keep your whole project in projected space.

ArcMap Projection double-check

I recently posted about how I had written a python script to convert Smith/Sandwell topography data to an ARC ASCII grid. I brought my resulting data in Arc and everything appeared to be correct. It all seemed to be going quite swimmingly; until that is, I noticed that the Prime Meridian appeared to pass just off the eastern coast of Australia. Hmmmm.. Last I checked, the Prime Meridian was still in Greenwich, UK. Clearly something was wrong. In my grid image, Greenwich fell on the very edge of my map. It seems that while I can bring the data into Arc and tell it what the projection is, Arc has a hard time displaying correct Lat/Long coordinates if you range of from 0 ->360 rather than -180 -> 180. There does not seem to be a way to specify that your longitude range is 0 -> 360 versus -180 -> 180 in Arc. I find this strange and am wondering if I am simply just missing it. Anyway, since I could not figure how to make it right in Arc, I went back to my original python code.

I decided to make my grid actually go from -180 -> 180 by amending my code. I basically grabbed the western half of my grid and appended it to the beginning of my file. I then cut off the redundant data. Using Struct.unpack to decode binary data results in a tuple, which cannot be modified, so the first step is to convert to a

#now unpack img data write out the ASCII file data
for j in range(nrow):
if v:
if j% 500 == 0:
print 'row:', j #print row # every 500 rows
raw_data = img.read(2*ncol)
row_tuple = struct.unpack('>'+str(ncol)+'h', raw_data)

#now move western half of data to the east
row = list(row_tuple)
east = row[ncol/2:ncol]
row.reverse()
east.reverse()
row.extend(east)
row.reverse()
corrected_row = row[0:ncol]

I had to reverse the rows before I could append the data, because the extend command places the data to be appended at the end of the row only. Once the data is correctly added, I simply reverse it back and snip off the extra. I could also just as easily grabbed the first half of the row instead (east = row[0:ncol/2]) and avoided the reversals, but this just happened to be how I worked through it first.

UPDATE: I showed Kurt my code last night and pointed out how I went about switching up the columns in my row variable. He said it was good that I figured out the two long ways, and then he showed me the one line method. Sure, you cannot really modify tuples after they have been created (e.g. no tuple.append or tuple.extend) but you can grab sections of them and switch them around like so:

In [1]: t = (1.2,3.4,5.6,7.8,9.1,10.2)

In [2]: t = t[3:] + t[:3]

In [3]: t
Out[3]:
(7.7999999999999998,
9.0999999999999996,
10.199999999999999,
1.2,
3.3999999999999999,
5.5999999999999996)

Applying this method to my code, I can switch up the western and eastern halves of the line I read in from the img file simply by:

#now move western half of data to the east
row = row[ncol/2:] + row[:ncol/2]

Pylint: The Python Code bug/quality checker

Now that I am off and running in Python, I thought I should install Pylint. Pylint will review your code for bugs and quality (as determined by some formalized Python standards).

You can run it via command line simply by typing:

pylint myscript.py

There are also a lot of options for silencing reports, getting more detail from messages, etc. If you leave the reports on, you get a global evaluation score. I just ran my script the decodes a binary SIO IMG file into an ARC ASCII grid through pylint and got this result:

Global evaluation
-----------------
Your code has been rated at 4.92/10

Not bad, considering this is my first ever script. My code runs fine; the majority of the messages involved formatting issues, such as some of my lines being too long. For code readability, the convention is to keep lines as short as possible. I definitely have some that I could break up into 2 lines. I also get some messages about not following my commas with a space (picky, picky).

After shortening just a couple lines and putting in some spaces after some commas, I ran it again:

Global evaluation
-----------------
Your code has been rated at 6.44/10 (previous run: 4.92/10)

It is still complaining about some line lengths and the fact that I have some variables named with single letters, but in general my code is pretty good. If you want to help make your code better, especially if you are just learning like I am or if you plan to distribute your code, I recommend Pylint.

Leaving the shallows...

So I recently have experienced a huge change in my research direction, and for those possibly going through the same thing, I thought I would share the hows and whys.

Last year, my original research project began taking a serious nose-dive. I knew what I wanted to do and thought I had a good plan of attack, but I simply could not account for some variables. Amongst other issues, there was some crucial information I was going to need to make it work, and I was not going to be able to get it. Lesson learned: beware projects involving potential proprietary information. However, if someone wants to give me hundreds of thousands of dollars so I can buy my own plane, lidar system, and hire a pilot, I may still be able to do it though...

Anyway, after a couple months of feeling like I was drowning in trying to find another project, I was getting scared. I was nearing the start of my third year and I had no project. Sure, I had like 45 credits worth of class that I had taken, but my research had basically gone up in smoke. I got to the point were I was considering taking a leave from school to give myself time to decide what I wanted to do and if I even wanted to continue in grad school. The best advice I got was to talk to some of the faculty outside my department, particularly young, fresh-faced, faculty that were still full of excitement. I did so, and it saved me.

Now, I have finally came across a project that I think is a perfect fit. It brings me back to my undergrad days of awe when I was first learning about oceanography and geology and everything amazed me. I am leaving the shallows and heading into the deep! My new research focuses (in a nutshell) on using high-resolution multibeam sonar data to determine fault structure and earthquake potential along oceanic transform faults. I am so excited. My committee is also very excited and very supportive. This is key, as last year this was not always the case.

I am sharing this because I think it is important for other grad students (especially the PhD-ers out there) to know that sometimes setbacks happen, and sometimes they are for the best. This new project definitely has me more excited than my old. There is much more room for collaboration, and I am going to get to work with some great folks down at Woods Hole! My committee changed slightly as a result of the new direction, and is more excited and supportive than ever before. Did I lose some time? Probably about 6 months of research time. Does it matter? No. In the long run, who cares? As one of my committee members said to me, "the setback itself matters not, what does matter is whether or not you bounce back." Setbacks in grad school are common. They happen to the best of us. After all, we are still learning how to ask the questions that lead to great projects in the first place.

So 2010 certainly has been a fresh start for me. A new year, a new project, and now, a new blog layout.

Converting SIO img file format to Arc ASCII grid w/ Python

Recently I have been battling with getting a SIO binary file format of predicted bathymetery (from Smith/Sandwell satellite topography data) into an Arc-friendly ASCII grid format. GMT has a nifty img2grd command which generates a netcdf of the img file; however, it takes some manual tweaking afterwards to get the coordinate bounds of the data to be Arc-happy (even if one follows img2grd grd2xyz with the -E option). I have decided I want a one-stop solution, where I could just feed in the img file and spit out an ARC ASCII grid and this means I need to write my own script. Therefore, I sat myself down this afternoon and began to teach myself Python (with some input from Kurt). Now, after just a couple hours of Googling, tweaking, and testing, I have a working Python script that does just what I need it to. It reads in an img file, and spits out a space-delimited ASCII grid file complete with the ARC header.

In order to test that I was decoding the binary properly, I wrote the first row of data out to its own little text file and used gnuplot to graph it up.

in my script I have the following:

for i in range(1):
row = struct.unpack('>'+str(ncol)+'h',img.read(2*ncol))
for depth in row:
    o.write(str(depth)+"\n")
print "\n"

then in terminal I call gnuplot and at the prompt type:

plot 'filename' with l

Certainly looks like a nice depth profile to me:

Org-mode for note organization

New Year: new attempt at note organization

I am always trying to find the best way to organize my random work notes and tips that I tend to collect in various notebooks, post-its, etc. I started keeping a text file last year, which was rather plain-Jane, but did the trick. Kurt, however, just introduced me to org-mode in Emacs, and I love it! Org-mode is extremely powerful. You can insert date stamps, create tags and link across multiple entries, create levels and sublevels and fold/unfold sections, highlight code blocks according to the scripting language, etc. Included html links are active from the get-go, and you can even include equations, graphics, and tables. One of the best parts is that you can export to HTML, PDF (via LaTeX), and DocBook.

Here is a sample of what my notes look like in Aquamacs using org-mode. I have organized my notes by program/tool, but you could also organize by date and use org-mode to keep a daily log. On the org-mode website, you'll see tons of examples. You can also use org-mode as a task scheduler and day planner. The possibilities seem endless.


And here is the resulting HTML file (ignore the old dates, I imported my old notes from last year):