get_coord_from_GeoTIFF

# -*- coding: utf-8 -*-
from osgeo import osr, gdal
import sys
import shapefile as sf
import math
import psycopg2

# 2013/09/03

# USAGE:
# python get_coord.py ps_RS2_20130829_072059_0076_SCWA_HHHV_SGF_279190_9158_8821989.zip.tif 20130829-20130830.txt

def create_table():
    try:
        conn = psycopg2.connect(database='gis', user='gisuser')
        curs = conn.cursor()
        # Try create db for IAPB data
        curs.execute("""CREATE TABLE my_drift (id integer, date1 date, date2 date, geom1 geometry, geom2 geometry)""")
        conn.commit()
        conn.close()
    except:
        print 'The table already exist or could not be created!'

def insert_data(iid,idate1,idate2,ilon1,ilat1,ilon2,ilat2):
    ''' Insert data into DB '''
    try:
        conn = psycopg2.connect(database='gis', user='gisuser')
        curs = conn.cursor()
        sql = "INSERT INTO my_drift (id,date1,date2,geom1,geom2) VALUES (%s,\'%s\',\'%s\',ST_GeomFromText(%s),ST_GeomFromText(%s));" % \
        (iid,idate1,idate2,'\'POINT( %s %s )\'' % (ilon1,ilat1),'\'POINT( %s %s )\'' % (ilon2,ilat2))
        #params = [int(id),dt_str,'POINT( %s %s )' % (lon,lat)]
        #print sql
        curs.execute(sql)
        conn.commit()
        conn.close()
        print 'OK! Data inserted!'
    except:
        print 'Insert error!'

def big_circle(llat1,llong1,llat2,llong2):
    #pi - число pi, rad - радиус сферы (Земли)
    rad = 6372795

    #в радианах
    lat1 = llat1*math.pi/180.
    lat2 = llat2*math.pi/180.
    long1 = llong1*math.pi/180.
    long2 = llong2*math.pi/180.

    #косинусы и синусы широт и разницы долгот
    cl1 = math.cos(lat1)
    cl2 = math.cos(lat2)
    sl1 = math.sin(lat1)
    sl2 = math.sin(lat2)
    delta = long2 - long1
    cdelta = math.cos(delta)
    sdelta = math.sin(delta)

    #вычисления длины большого круга
    y = math.sqrt(math.pow(cl2*sdelta,2)+math.pow(cl1*sl2-sl1*cl2*cdelta,2))
    x = sl1*sl2+cl1*cl2*cdelta
    ad = math.atan2(y,x)
    dist = (ad*rad) #/1000. [m]

    #вычисление начального азимута
    x = (cl1*sl2) - (sl1*cl2*cdelta)
    y = sdelta*cl2
    z = math.degrees(math.atan(-y/x))

    if (x < 0):
        z = z+180.

    z2 = (z+180.) % 360. - 180.
    z2 = - math.radians(z2)
    anglerad2 = z2 - ((2*math.pi)*math.floor((z2/(2*math.pi))) )
    angledeg = (anglerad2*180.)/math.pi

    print 'Distance >> %.2f' % dist, ' [km]'
    print 'Initial bearing >> %.0f ' % angledeg, '[degrees]\n'
    return '%.0f' % dist, '%.0f' % angledeg

def read_res_file(res_file):
    ''' Readlines of out file with vectors '''
    ff = open(res_file)
    lines = ff.readlines()
    # Exclude 1st str
    lines = lines[1:]

    # Start make shape here
    shp_name = res_file.split('.')[0]
    shp_filename = "shp/" + shp_name
    #create the shapefile
    w = sf.Writer(sf.POLYLINE)
    w.field('vector_num','C','40')
    w.field('Magnitude','C','40')
    w.field('Azimuth','C','40')
    i = 0
    # Create txt file
    ff_txt = open('txt/%s.txt' % shp_name,'w')
    for line in lines:
        xx1 = line.split()[0]
        yy1 = line.split()[1]
        xx2 = line.split()[2]
        yy2 = line.split()[3]

        print 'x1: %s' % xx1
        print 'y1: %s' % yy1
        x1 = gt[0] + float(xx1)*pixelWidth
        y1 = gt[3] + float(yy1)*pixelHeight
        latlon = transform.TransformPoint(float(x1),float(y1))
        print latlon
        lon1 = latlon[0]
        lat1 = latlon[1]

        print 'x2: %s' % xx2
        print 'y2: %s' % yy2
        x2 = gt[0] + float(xx2)*pixelWidth
        y2 = gt[3] + float(yy2)*pixelHeight
        latlon = transform.TransformPoint(float(x2),float(y2))
        print latlon
        # Length in [meters]
        length = math.hypot(x2-x1,y2-y1)
        print 'Length: %.2f [km]' % (math.hypot(x2-x1,y2-y1)/1000.)
        lon2 = latlon[0]
        lat2 = latlon[1]

        w.line(parts=[[[lon1, lat1],[lon2, lat2]]])

        # Big circle length
        mag, az = big_circle(float(lat1),float(lon1),float(lat2),float(lon2))

        w.record(str(i),str(mag),str(az))

        ff_txt.write('%5.2f %5.2f %5.2f %5.2f\n' % (lon1,lat1,lon2,lat2))

        # 2DB
        idate1 = res_file.split('.')[0][0:4] + '-' + res_file.split('.')[0][4:6] + '-' + res_file.split('.')[0][6:8]
        idate2 = (res_file.split('.')[0]).split('-')[-1][0:4] + '-' + (res_file.split('.')[0][4:6]).split('-')[-1] + '-' + (res_file.split('.')[0]).split('-')[-1][6:8]
        insert_data(i,idate1,idate2,lon1,lat1,lon2,lat2)

        i = i + 1
    ff.close()
    ff_txt.close()

    # Save shape
    w.save(shp_filename)
    # create the PRJ file
    prjfn = shp_filename[:]
    prj = open("%s.prj" % prjfn, "w")
    epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
    prj.write(epsg)
    prj.close()

    return latlon

# get the existing coordinate system
in_file = sys.argv[1]
res_file = sys.argv[2]

ds = gdal.Open(in_file)
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs.ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs)

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()

originX = gt[0]
originY = gt[3]
pixelWidth = gt[1]
pixelHeight = gt[5]


minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5]

#get the coordinates in lat long

latlon=read_res_file(res_file)

#latlong = transform.TransformPoint(3725,32)