Showing posts with label SIO img. Show all posts
Showing posts with label SIO img. Show all posts

Thursday, January 21, 2010

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]

Sunday, January 17, 2010

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: