Cataclysmic Mutation

Machine Learning and Whatever Else

Geotagging in C (Part 2: NMEA)

When last we left off, we had written most of the code necessary for reading and writing the metadata to and from Nikon raw images in NEF format. Recall, the ultimate goal is to be able to match the timestamp from a photo to an entry from our GPS log and write the corresponding location data back into the EXIF data for the photo. In this installment, we’ll look at how to parse the needed information from the files output by the data logger.

If you recall, the GPS unit I bought (the AMOD AGL3080), recorded its log files in a format called NMEA. NMEA takes its name from the National Marine Electronics Association, and it is, in a word, bizarre. Note that NMEA in actuality is not simply a data format, but defines an entire communications protocol. However, since the simple data logger we’re interested in isn’t really communicating with anything, all we really need to know is how to parse a file of data that corresponds to the NMEA protocol.

.NMEA data is at heart nothing more than CSV formatted data, but in NMEA, there are different types of “sentences”. Each line of data begins with a token that defines the sentence type, and the sentence type tells you what the remaining fields in the CSV formatted line represent. A reasonably complete listing of NMEA sentences and their meaning can be found here. A not so quick glance through that information shows that we need at least the NMEA sentence type “$GPGGA”. The GPGGA format is reproduced below.

$GPGGA,hhmmss.ss,llll.ll,a,yyyyy.yy,a,x,xx,x.x,x.x,M,x.x,M,x.x,xxxx*hh
1    = UTC of Position
2    = Latitude
3    = N or S
4    = Longitude
5    = E or W
6    = GPS quality indicator (0=invalid; 1=GPS fix; 2=Diff. GPS fix)
7    = Number of satellites in use [not those in view]
8    = Horizontal dilution of position
9    = Antenna altitude above/below mean sea level (geoid)
10   = Meters  (Antenna height unit)
11   = Geoidal separation (Diff. between WGS-84 earth ellipsoid and
       mean sea level.  -=geoid is below WGS-84 ellipsoid)
12   = Meters  (Units of geoidal separation)
13   = Age in seconds since last update from diff. reference station
14   = Diff. reference station ID#
15   = Checksum

At first glance, it would appear that this is all the information we need. It contains a timestamp, the latitude, longitude, and some information that tells us how reliable the GPS information was for this particular fix. However, there is a significant problem, which will be evident when we see an example of a real GPGGA sentence.

$GPGGA,111929.943,6400.0478,N,02058.8074,W,1,06,1.6,44.6,M,60.4,M,,0000*7B

The record show above tells us that at 11:19:29.943 UTC, the unit was at 6400.0478 north latitude, and 02058.8074 west longitude. Ignoring for a moment that the latitude and longitude seem to make no sense (we’ll come back to that in a moment), the bigger problem is that there is no date information. We know the time, but not the date. If we had taken photos on more than one day, there would be no way to figure out which GPGGA record we should match with each photo.

Looking back at the NMEA spec, it looks like the “$GPRMC” sentence contains the date information. Note that it also contains the latitude, longitude, and time information again, so if we’re content with that, we could simply dispense with the GPGGA sentences altogether. However, GPGGA sentences contain the secondary GPS fix information, as well as the altitude, and I wanted to be able to record altitude if I knew it. Thus, the strategy seems to be this: for each picture, we’ll find both the GPGGA and GPRMS records nearest to the time the photo was taken and combine their data to generate the EXIF data to write into the photo. Here is where I went ahead and cheated. The GPS unit I bought always seems to record GPGGA and GPRMS records every second, so I simply read each GPGGA sentence, left the date blank, and then filled it in based on the next GPRMS sentence read from the file. I can’t guarantee this is a bug-free strategy in general, but it does appear to work for me.

Now let’s return to the bizarrely formatted latitude and longitude for a moment. If you think back to your high school geometry class, there are several ways to represent portions of a circle. With latitude and longitude, the conventional units are degrees measured from specified zero points (the equator being zero degrees latitude, and zero degrees longitude being set through Greenwich). Again, there are a couple of ways to specify coordinates: we can specify a triplet of degrees, minutes, and seconds, or we can specify decimal degrees. For example, two latitudes 46º30’15” is the same point as 46.504167º. We get this by realizing that there are 60 minutes in a degree and 60 seconds in a minute, so 46º30’15” becomes 46 (30/60) (15/3600)=46.504167º.

So which format does NMEA choose? Both…neither…well, it’s complicated. NMEA, or at least as implemented by my device, chooses to store coordinates as degrees and decimal minutes. Oh, and the decimal point is in the wrong place. This is presumably what one gets when one lets fishermen design a data storage and transfer format. In any case, a coordinate like 6400.0478 N, 02058.8074 W from our example above actually represents 64º00.0478’ north latitude, 20º58.8074’ west longitude. Once we know this, it’s easy enough to convert to the degrees-minutes-seconds format needed for the EXIF convention. Note that EXIF also does not allow the character mnemonics like ‘N’ and ‘W’, so we’ll need to detect south and east coordinates and express them as negative numbers.

Let’s get some of this down in code. First, we need a structure to keep track of the notion of a timestamped location. Recall, we’ll fill this from both a GPGGA and a GPRMS sentence from the NMEA data log.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
typedef struct {
    time_t when;
    char status; /* A (active), V (void) */
    double latitude;
    char lat_ref;
    double longitude;
    char lon_ref;
    double speed; /* in knots */
    double heading; /* in degrees */
    double altitude; /* in meters */
    double geoid_ht; /* in meters */
    int num_sat;
    int quality;
} location_t;

We had some choices to make here. Do we store the directional mnemonics or convert the degrees as we read the file? Do we store the information that doesn’t seem to be immediately useful? What format should we store the time in? Some of these questions are somewhat arbitrary. The choices I made weren’t strongly considered; they just seemed reasonable to me. I store the time and date information in a standard time_t value, primarily because we can easily compare values allowing an efficient binary search of our array of locations read from the file. I decided to store some information from the NMEA log that didn’t seem necessary (heading, number of satellites, etc.). This was for no real reason other than interest – it was data that I thought I might want to do something with.

As NMEA files are CSV files, we need to be able to parse those files. Again, as with TIFF files and NMEA files, there are libraries available to handle this for us, but the whole point of this project was to write some fun low-level code, so I chose to implement my own simple CSV file parser.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/*
 * parse a line of comma separated text into a given
 * array of tokens
 */
void parse_line(char* line, char* sep, char** toks)
{
    char*  curr;
    char*  p;
    int i = 0;

    i=0;
    while(line != NULL && i < NUM_TOKENS) {
        curr = strsep(&line, sep);

        /* NMEA files from my device are DOS formatted; kill the rn stuff */
        if((p = strchr(curr, 'r')) != NULL)
            *p = '';
        else if((p = strchr(curr, 'n')) != NULL)
            *p = '';

        /* note we're reusing the space in the line here */
        toks[i++] = curr;
    }
}

With that code, we can read an entire NMEA file as follows.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
/*
 * parse a file of NMEA sentences into an array of location_t records
 */
void parse_nmea_file(FILE* fp, location_t** rows, int* num_recs, int max_size)
{
    char** toks;
    char*  line;

    /* init some memory to use to parse the nmea file */
    toks = (char**)malloc(NUM_TOKENS * sizeof(char*));

    /* parse each GPRMC record */
    *num_recs = 0;
    line = (char*)malloc(82 * sizeof(char));
    while((line = fgets(line, 82, fp)) != NULL) {
        parse_line(line, ",", toks);
        if(strncmp(toks[0], "$GPRMC", MAX_TOKEN_LEN) == 0) {
            init_rmc_rec(&(*rows)[(*num_recs)  ], toks);
        } else if(strncmp(toks[0], "$GPGGA", MAX_TOKEN_LEN) == 0) {
            /* if first record is a GPGGA record, ignore it */
            if(*num_recs > 0)
                process_gga_rec(&(*rows)[(*num_recs)-1], toks);
        } else {
            /* skip other sentence types */
            continue;
        }

        /* if we've run out of space, realloc twice the space and keep going */
        if(*num_recs >= max_size) {
            max_size = max_size * 2;
            *rows = (location_t*)realloc(*rows, max_size * sizeof(location_t));
            if(!*rows) {
                fprintf(stderr, "error enlarging nmea sentences arrayn");
                exit(EXIT_FAILURE);
            }
        }
    }
    free(toks);
    free(line);
}

The strategy here is quite simple. We allocate an initial array to hold the location records and begin reading lines. When we fill up the array, we realloc some more memory and continue reading. This is a classic strategy for implementing growable arrays. As each line is read, we check the first token in the sentence. If it’s a GPRMC sentence, we initialize a new location entry. We parse the time and date and pull out the GPS information that is present. The rest we set to zeros. As soon as we read a GPGGA record, we pass the same partially populated location record in and fill in the remainder of the GPS fix information. All other sentence types are ignored. Note that the latitude and longitude are stored in the NMEA native degrees/decimal minutes format.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
/* 
 * builds up a location_t record based on a parsed GPRMC sentence
 */
void init_rmc_rec(location_t* rec, char** toks)
{
    struct tm fix_time;
    char tmp[3];

    /* time */
    strncpy(tmp, toks[1], 2);
    fix_time.tm_hour = atoi(tmp);
    strncpy(tmp, toks[1] 2, 2);
    fix_time.tm_min = atoi(tmp);
    strncpy(tmp, toks[1] 4, 2);
    fix_time.tm_sec = atoi(tmp);

    /* date */
    strncpy(tmp, toks[9], 2);
    fix_time.tm_mday = atoi(tmp);
    strncpy(tmp, toks[9] 2, 2);
    fix_time.tm_mon = atoi(tmp) - 1;
    strncpy(tmp, toks[9] 4, 2);
    fix_time.tm_year = 2000   atoi(tmp) - 1900;

    /* set the offset from UTC to 0, as GPS reports times in UTC anyway */
    fix_time.tm_gmtoff = 0;

    rec->when = timegm(&fix_time);
    rec->status = toks[2][0];
    rec->latitude = atof(toks[3]);
    rec->lat_ref = toks[4][0];
    rec->longitude = atof(toks[5]);
    rec->lon_ref = toks[6][0];
    rec->speed = atof(toks[7]);
    rec->heading = atof(toks[8]);

    /* these fields aren't part of GPRMS sentence; we'll add them later */
    rec->altitude = 0;
    rec->geoid_ht = 0;
    rec->quality = 0;
    rec->num_sat = 0;
}

/*
 * updates a given location_t record with information from a parsed GPGGA sentence
 */
void process_gga_rec(location_t* rec, char** toks)
{
    /* if we've already seen a more accurate GPGGA record and used it,
       then bail.  this is because there may be multiple GPGGA records
       per GPRMC record, and we want to use the first one succeeding a
       GPRMC record, as it should have the closest timestamp. */
    if(rec->altitude > 1e-3)
        return;

    /* make sure we have at least some kind of fix */
    if(atoi(toks[6]) == 0)
        return;

    /* update the previous location_t record with the fix and altitude data
       from the currently parsed GPGGA sentence tokens */
    rec->quality = atoi(toks[6]);
    rec->num_sat = atoi(toks[7]);
    rec->altitude = atof(toks[9]);
    rec->geoid_ht = atof(toks[11]);
}

Once this process is complete and the file has been completely read, what we’re left with is an array of location_t structures, ordered by time. Note that if we wanted to be completely general, we should go ahead and sort the array, but if we only read one GPS log file during the run, it will be in order by time already. Since the file is ordered, we can implement a simple binary search based on time.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
/*
 * custom binary search that looks for the record in the array whose timestamp
 * is nearest to the given timestamp.  if no records are stamped within epsilon
 * seconds, returns NULL.
 */
location_t* find_location_at(location_t* rows, unsigned int nrows, time_t ts, int epsilon)
{
    int low = 0;
    int mid;
    int high = nrows - 1;

    while(low <= high)
    {
        mid = (low + high) / 2;
        if(rows[mid].when < ts)
            low = mid + 1;
        else if(rows[mid].when > ts)
            high = mid - 1;
        else
            return &rows[mid];
    }

    /* low has passed high, so check which is closer to desired time */
    if(fabs(rows[high].when - ts) <= fabs(rows[low].when - ts) &&
       fabs(rows[high].when - ts) <= epsilon)
        return &rows[high];
    else if(fabs(rows[low].when - ts) < fabs(rows[high].when - ts) &&
       fabs(rows[low].when - ts) <= epsilon)
        return &rows[low];

    /* no record found within required time limit */
    return NULL;
}

Note that there’s really no excuse, even for a project where the point was to have fun writing low-level code, to avoid the built-in bsearch function from the C library. However, we want slightly different behavior from a standard binary search. In particular, it’s possible that when you’re out walking around taking photos, your GPS fix may drop out for a few minutes. Tall buildings can cause problems in cities, or maybe you simply forgot to turn it on for a few minutes when you started your photo walk. The nature of geotagging is such that we do not demand perfection. If the choice is between not tagging a photo or tagging it incorrectly but within a few hundred feet, most people would choose the approximately correct location. In any case, what we want is the ability to specify a tolerance – the largest time difference that we’re willing to accept as matching a given photo. A standard binary search will find an element if it exists, but has no way to say, “find the closest element unless the distance between it and the probe is greater than some value”. The custom binary search shown above allows just that.

So we can now read an NMEA log, store an array of locations with timestamps, and search this array for the best available location given a probe timestamp. From the previous post, we have the ability to parse information from an NEF file. What’s left is the specific code to get the capture time from a set of images, to wire up the search to find the best location, and to write the location back into the photo. We’ll pick up there in the next installment.