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.

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:

Perls, Pythons, and Rs...

I am not a coder, plain and simple. I may work in Matlab some, but even that involves a lot of trial and error for me. What would take some of my friends 5 minutes can take me an hour or more as I muddle my way through function calls and syntax. I am definitely keen to learn though. In my geodesy class we had a semester-long project in which we became a virtual GPS receiver and had to solve for different position fixes (basically using multiple iterations of a least-squares analysis). Most people pick Excel to do it in, including the prof. All that copy and paste and only one line in which to type your equations just gets me, and I cannot do it. I decided on Matlab. I definitely struggled with it a lot, but it was relatively easy to debug when there were issues, and in the end it actually worked.

Now one of my other profs is trying to sway me towards R instead of Matlab. The draw of R is that it is freeware, highly supported by users, and completely object-oriented. Apparently it is rapidly becoming the main coding software of statisticians. The drawback, for me obviously, is switching over and starting from scratch all over again. Other profs tell me to stick with Matlab, it is more mainstream and therefore more widely-accepted. I wonder though, which is better? Perhaps R is more intuitive, which would certainly be nice...

Perl versus Python is another deciding point for me. I recently worked with some Python code involving the same tide data from the blog below. Python seems relatively straightforward on first glance and is pretty powerful. It certainly ran through an entire month of 12-second tide data and converted all the raw temperatures, water levels, and time stamps in mere seconds. Mightily impressive. Then someone told me about Perl. It is also a great language and can do the same things as Python. If you look it up online, people rant and rave about both. A lot of folks seem to think that while Python is somewhat easier to script in than Perl, Perl is still the better choice. Python appears to be more user-friendly, but Perl can work with other programming languages (such as C++) better than Python. Also, it seems that while Perl can do everything Python can, the vice-versa is not necessarily true. My question is for someone who really has no experience with either, but wants to learn one to help with data processing, which one is better? Should I go for the one that is easier to use, or the one that has more compatibility? Does it even really matter for what I will be doing? How much coding will analyzing full-waveform lidar data even require?

Blah, why are there so many options?