|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +import sys |
| 4 | +from os.path import isfile, sep |
| 5 | +from numpy import double, floor |
| 6 | + |
| 7 | +def translate(lat, lon, latdelta, londelta): |
| 8 | + latidx = int(floor((90. - lat) / latdelta + 1)) |
| 9 | + lonidx = int(floor((lon + 180.) / londelta + 1)) |
| 10 | + |
| 11 | + lat1 = 90. - latdelta * (latidx - 1) |
| 12 | + lat2 = 90. - latdelta * latidx |
| 13 | + lon1 = -180. + londelta * (lonidx - 1) |
| 14 | + lon2 = -180. + londelta * lonidx |
| 15 | + |
| 16 | + if lat > lat1 or lat < lat2 or lon < lon1 or lon > lon2: |
| 17 | + raise Exception('Not in cell!') |
| 18 | + |
| 19 | + return latidx, lonidx |
| 20 | + |
| 21 | +if len(sys.argv) != 6: |
| 22 | + print 'Usage: trans_gidx.py latidx1 lonidx1 lat_delta1[,lon_delta1] lat_delta2[,lon_delta2] weathdir2' |
| 23 | + sys.exit(1) |
| 24 | + |
| 25 | +latidx = int(sys.argv[1]) |
| 26 | +lonidx = int(sys.argv[2]) |
| 27 | +delta1 = [double(d) / 60. for d in sys.argv[3].split(',')] # arcminutes -> degrees |
| 28 | +delta2 = [double(d) / 60. for d in sys.argv[4].split(',')] |
| 29 | +weathdir2 = sys.argv[5] |
| 30 | + |
| 31 | +if len(delta1) == 1: |
| 32 | + latdelta1 = londelta1 = delta1[0] |
| 33 | +else: |
| 34 | + latdelta1, londelta1 = delta1 |
| 35 | +if len(delta2) == 1: |
| 36 | + latdelta2 = londelta2 = delta2[0] |
| 37 | +else: |
| 38 | + latdelta2, londelta2 = delta2 |
| 39 | + |
| 40 | +latrat, lonrat = latdelta2 / latdelta1, londelta2 / londelta1 |
| 41 | + |
| 42 | +latd = [0.5, 0, 1, 0, 1] # start at center |
| 43 | +lond = [0.5, 0, 0, 1, 1] |
| 44 | +latd += [0.5 + i for i in [0, 0, latrat, -latrat, latrat, latrat, -latrat, -latrat]] |
| 45 | +lond += [0.5 + i for i in [lonrat, -lonrat, 0, 0, lonrat, -lonrat, lonrat, -lonrat]] |
| 46 | + |
| 47 | +lats = [90. - latdelta1 * (latidx - i) for i in latd] |
| 48 | +lons = [-180. + londelta1 * (lonidx - i) for i in lond] |
| 49 | + |
| 50 | +for i in range(len(lats)): |
| 51 | + latidx2, lonidx2 = translate(lats[i], lons[i], latdelta2, londelta2) |
| 52 | + latidx2, lonidx2 = '%03d' % latidx2, '%03d' % lonidx2 |
| 53 | + wfile = sep.join([weathdir2, latidx2, lonidx2, '%s_%s.psims.nc' % (latidx2, lonidx2)]) |
| 54 | + if isfile(wfile): |
| 55 | + print latidx2 + sep + lonidx2 |
| 56 | + sys.exit(0) # success |
| 57 | + |
| 58 | +sys.exit(1) # failure |
0 commit comments