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]

1 comment: