GenericMappingTools/gmt

Calculate gradient and azimuth of a profile.

Esteban82 opened this issue · 5 comments

Esteban82 commented

For the work I am going to present at the AGU I made a small software (in C) that from a profile (with XYZ data) calculates the intermediate position between two successive data (Xm, Ym), and gradient and azimuth (Xm, Ym, Az, Grad). This way I can use it as input for greenspline (-A+f2). I think it would be nice to add it as a feature for GMT.

What do you think? My idea is to do it but I probably need help.

I think it should be added as a new argument for mapproject? Do you agree? What letter should I use?

This is the code:

#include <math.h>
#include <stdio.h>

int main() {
    // Definir los nombres de los archivos de entrada y salida
    const char* trackFileName = "Bad_Track";
    const char* outFileName = "Track_Gradient.txt";

    // Abrir los archivos de entrada y salida
    FILE *trackFile = fopen(trackFileName, "r");
    FILE *outFile = fopen(outFileName, "w");

    if (trackFile == NULL || outFile == NULL) {
        perror("Error al abrir archivos");
        return 1;
    }

    // Variables para almacenar los puntos medios
    double Xm, Ym, DX, DY, DZ, Distance, Gradient, Azimuth, Radianes;

    // Variables para leer las coordenadas de cada línea
    double x1, y1, z1, x2, y2, z2;

    // Leer las coordenadas de la primera línea
    if (fscanf(trackFile, "%lf %lf %lf", &x1, &y1, &z1 ) != 3) {
        perror("Error al leer el archivo de entrada");
        return 1;
    }

    // Loop para procesar cada línea
    while (fscanf(trackFile, "%lf %lf %lf" , &x2, &y2, &z2) == 3) {
        // 1. Calcular el punto medio
        Xm = (x1 + x2) / 2.0;
        Ym = (y1 + y2) / 2.0;

        // 2. Gradiente
        DX= (x2-x1);
        DY= (y2-y1);
        Distance= sqrt(pow(DX, 2) + pow(DY, 2));
        DZ = (z2-z1);
        Gradient= (DZ / Distance);
       
        // 3. Azimuth
        Radianes = atan2 (DX, DY);
        Azimuth = fmod ((Radianes * (180/M_PI)+360),360);
        //Azimuth = fmod (Azimuth * (180 /  M_PI) + 180, 360);

        // Escribir el punto medio en el archivo de salida
        fprintf(outFile, "%.2f %.2f %.5f %.4f\n", Xm, Ym, Gradient, Azimuth);

        // Actualizar las coordenadas anteriores
        x1 = x2;
        y1 = y2;
        z1 = z2;
    }

    // Cerrar los archivos
    fclose(trackFile);
    fclose(outFile);

    printf("Proceso completado con éxito.\n");

    return 0;
}
joa-quim commented
Esteban82 commented

Hi Joaquim

I understand that your suggestion is to add a new modifier (-f6) for greenspline in which one enters XYZ data and GMT internally converts that data to the -f2 format.

joa-quim commented

Whithout recheking greenspline docs, I would say, yes.

Esteban82 commented

Search also for functions that compute the azimuth instead of reimplement it.

Just to be sure. I should look up where the formulas are and call them directly. Not copy them, right?

For example, I understand that the azimuth is calculated inside mapproject (in line 504 I think). So I need to calculate that formula.

Then, I would need to do the same to calculate distances and gradient.

joa-quim commented

Yes, find the function. A not obvious exercise. Often I need to do debugging to find what I'm looking for.

This one may be a possible target but code in GMT do not call it directly. Instead this is called via a pointer to function where selection was done earlier based on the geog/geodesic model selected.

gmt/src/gmt_map.c

Line 2466 in da7f9a8

GMT_LOCAL double gmtmap_az_backaz_cartesian (struct GMT_CTRL *GMT, double lonE, double latE, double lonS, double latS, bool baz) {