Calculate gradient and azimuth of a profile.

[Esteban82]

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]

Olá Federico

In between flights so I’ll kill some time trying to reply to this from the mobile.

First, good to see your efforts diving into C.

Now, to the particular subject. I don’t think it would be a good idea to add that to mapproject that is a module already overloaded. Normally I would not be favorable to add this type of calculations as a per se but since you want to use it in greenspline, it is there, under the -A option, that it could be added.

On a more general comment/revision, reading from files (and output) is more complex than that. We must use the C lib IO functions to do that (that us, no direct fopen(...), fprintf, etc...). That is how we unsure that GMT continues to behave as a library.

Since the -A option already ingests a file, must see how to intercept that to compute the new variables and make space for them. Search also for functions that compute the azimuth instead of reimplement it.

Good luck

Joaquim

Enviado a partir do Outlook para iOShttps://aka.ms/o0ukef

---

De: Federico Esteban @.> Enviado: Tuesday, August 13, 2024 3:22:03 AM Para: GenericMappingTools/gmt @.> Cc: Subscribed @.***> Assunto: [GenericMappingTools/gmt] Calculate gradient and azimuth of a profile. (Issue #8565)

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

include

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;

}

— Reply to this email directly, view it on GitHubhttps://github.com/GenericMappingTools/gmt/issues/8565, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAEDF2LGX3PPGFRJ6M3AYE3ZREKOXAVCNFSM6AAAAABMM33HCWVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQ3DCOBSGM2DMMA. You are receiving this because you are subscribed to this thread.Message ID: @.***>

[Esteban82]

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]

Whithout recheking greenspline docs, I would say, yes.

[Esteban82]

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]

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. https://github.com/GenericMappingTools/gmt/blob/da7f9a8a28d46fef5cd9dc8c4ec10ebfddbdbf0b/src/gmt_map.c#L2466