oceanhackweek / ohw23_proj_drone_georef

We aim to develop a workflow to match drone images with a satellite grid
Apache License 2.0
1 stars 0 forks source link

Extract image positions into Open done map format #4

Open NickMortimer opened 1 year ago

NickMortimer commented 1 year ago

I have some code to do this https://docs.opendronemap.org/geo/

NickMortimer commented 1 year ago
""" 
    Column 1: Image sequence.
    Column 2: Second of week of the exposure timestamp in UTC time, with seconds expressed in GPS
    time format.
    Column 3: GPS week of the exposure timestamp in UTC time, with seconds expressed in GPS time
    format.
    Column 4: Offset (in mm) of the antenna phase center to camera CMOS sensor center in the north
    (N) direction at the time of exposure of each photo. It is positive when the CMOS center is in the
    north direction of the antenna phase center and negative when in the south direction.
    Column 5: Offset (in mm) of the antenna phase center to camera CMOS sensor center in the east
    (E) direction at the time of exposure of each photo. It is positive when the CMOS center is in the
    east direction of the antenna phase center and negative when in the west direction.
    Column 6: Offset (in mm) of the antenna phase center to camera CMOS sensor center in the vertical
    (V) direction at the time of exposure of each photo. It is positive when the CMOS center is below
    the antenna phase center and negative when the former is above the latter.
    Column 7: real-time position latitude (Lat) of the CMOS center acquired at the time of exposure,
    in degrees. When the aircraft is in the RTK mode, its position is the RTK antenna phase center
    position plus the offset of the antenna phase center relative to the CMOS center at the time of
    exposure, with the RTK accuracy (centimeter level); and when the aircraft is in the GPS mode,
    its position is that detected by GPS single-point positioning plus the offset of the RTK antenna
    phase center relative to the CMOS center at the time of exposure, with GPS single-point
    positioning accuracy (meter level).
    Column 8: real-time position longitude (Lon) of the CMOS center acquired at the time of exposure,
    in degrees.
    Column 9: real-time height of the CMOS center acquired at the time of exposure, in meters. The
    height is a geodetic height (commonly known as the ellipsoid height), relative to the surface of
    the default reference ellipsoid (WGS84 is the default and can be set as other ellipsoids such as
    CGCS2000 via a different CORS station system/reference). Note that the above height is not based
    on the national 85 elevation datum or 56 elevation datum (normal height) or commonly used
    EGM96/2008 elevation datum (orthometric height) worldwide.
    For correlation of orthometric height, normal height, and geodetic height, refer to the following link:
    http://www.esri.com/news/arcuser/0703/geoid1of3.html
    Columns 10 to 12:
    Standard deviation (in meters) of positioning results in the north, east, and the vertical direction,
    describe the relative accuracy of positioning in the three directions.
    Column 13:
    RTK status flag. 0: no positioning; 16: single-point positioning mode; 34: RTK floating solution;
    50: RTK fixed solution. When the flag bit of a photo is not 50, it is recommended not to use this
    photo directly in map building.
    Rinex.Obs file: real-time decoded satellite observation file (GPS+GLO+BDS+GAL) received by
    the UAV, in RINEX 3.02 format. This file can be directly imported into PPK post-processing
    software for post-processing.

    1   268619.354051   [2129]      24,N         9,E       193,V    -21.97273819,Lat    113.93433129,Lon    91.071,Ellh 1.056360, 1.085926, 3.245131    16,Q
    2   268625.283927   [2129]      20,N        15,E       193,V    -21.97253139,Lat    113.93438555,Lon    91.182,Ellh 1.097883, 1.051731, 3.346695    16,Q
    3   268627.779356   [2129]      10,N        13,E       194,V    -21.97234525,Lat    113.93440346,Lon    91.287,Ellh 1.095263, 1.051843, 3.344983    16,Q
    4   268630.314319   [2129]      13,N        10,E       194,V    -21.97216604,Lat    113.93443131,Lon    89.959,Ellh 1.090093, 1.041215, 3.321528    16,Q
    5   268632.830491   [2129]      15,N         6,E       194,V    -21.97199354,Lat    113.93445195,Lon    89.447,Ellh 1.081201, 1.031734, 3.327098    16,Q
    6   268635.351011   [2129]      14,N         8,E       194,V    -21.97182339,Lat    113.93448809,Lon    92.148,Ellh 1.076505, 1.022308, 3.273235    16,Q
    7   268638.543134   [2129]      66,N        42,E       178,V    -21.97168707,Lat    113.93449202,Lon    88.557,Ellh 1.110462, 1.036878, 3.367443    16,Q
    8   268642.289662   [2129]      20,N         3,E       194,V    -21.97167044,Lat    113.93429231,Lon    91.825,Ellh 1.109429, 1.020208, 3.272965    16,Q

    https://hpiers.obspm.fr/eoppc/bul/bulc/Leap_Second.dat
"""

import pandas as pd
import os
import glob
import ftplib
import re

"""
    Get 30 second data
"""
import time
secsInWeek = 604800 
secsInDay = 86400 
gpsEpoch = (1980, 1, 6, 0, 0, 0)  # (year, month, day, hh, mm, ss) 

def UTCFromGps(gpsWeek, SOW, leapSecs=16): 
    """converts gps week and seconds to UTC 

    see comments of inverse function! 

    SOW = seconds of week 
    gpsWeek is the full number (not modulo 1024) 
    """ 
    secFract = SOW % 1 
    epochTuple = gpsEpoch + (-1, -1, 0)  
    t0 = time.mktime(epochTuple) - time.timezone  #mktime is localtime, correct for UTC 
    tdiff = (gpsWeek * secsInWeek) + SOW - leapSecs 
    t = t0 + tdiff 
    #print t
    (year, month, day, hh, mm, ss, dayOfWeek, julianDay, daylightsaving) = time.gmtime(t) 
    #use gmtime since localtime does not allow to switch off daylighsavings correction!!! 
    return (year, month, day, hh, mm, ss + secFract) 

def read_mrk(filename,leapseconds=37,index='Sequence'):
    names=['Sequence','GPSSecondOfWeek','GPSWeekNumber','NorthOff','EastOff','VelOff','Latitude','Longitude','EllipsoideHight','Error','RTKFlag']
    data = pd.read_csv(filename,header=None,sep='\t',names=names)
    if data.shape[0]>0:
        data['Latitude']=data.Latitude.str.split(',',expand=True)[0].astype(float)
        data['Longitude']=data.Longitude.str.split(',',expand=True)[0].astype(float)
        data['GPSWeekNumber']=data.GPSWeekNumber.str[1:-1].astype(int)
        data['UTCTime'] =data.apply(lambda x: pd.to_datetime("1980-01-06 00:00:00")+pd.Timedelta(weeks=x['GPSWeekNumber'])+pd.Timedelta(seconds=x.GPSSecondOfWeek)
                                    -pd.Timedelta(seconds=leapseconds),axis=1)
        data.set_index(index,inplace=True)
    return data 
def read_mrk_gpst(filename,index='UTCTime'):
    return read_mrk(filename,leapseconds=0,index=index)

def read_pos(filename):
    # Using readlines()
    with open(filename, 'r') as pos:
        head = 0
        while pos:
            line = pos.readline()
            head= head+1
            if line.startswith(r'%  GPST'):
                break
    names = re.split('\s+',line.strip())
    data = pd.read_csv(filename,skiprows=head,sep='\s+',names=names,header=None)
    data.index =pd.to_datetime(data['%']+' '+data['GPST'])
    data.drop(['%','GPST'],axis=1,inplace=True)
    data.index.name = 'GPST'
    return data

def imagenames(filename):
    path =os.path.split(filename)[0]
    base =os.path.split(filename)[1].split('Timestamp')[0]+'*.JPG'
    images = glob.glob(os.path.join(path,base))
    data = pd.DataFrame({'Path':images})
    data['Sequence'] = data.Path.str.split('_').str[-1].str.split('.').str[0].astype(int)
    data.set_index('Sequence',inplace=True)
    return data
alsonathif commented 1 year ago

image positions added to open drone map format as 'gcp_list.txt'