Cześć, szukałem pomocy w różnych miejscach, pisałem do kilkunastu osób na olx, ale najczęstsza odpowiedź: C, GDAL? sorry, nie zajmuję się tym. Dlatego mam nadzieję, że tutaj znajdę pomoc :)
Mam kod który jest jaki jest, ale za to nie działa. Uruchamia się, woła o pliki wejściowe (.tif i współrzędne.txt), przetwarza i daje plik wynikowy. W teorii super, a w praktyce daje złe wyniki w pliku wynikowym. Dlaczego? Prawdopodobnie ma problem z Geo transform, ale tego nie jestem pewien, bo sam nie ogarniam C. Działam w PuTTY, Midnight Commander.
#include <stdio.h>
#include <stdlib.h>
#include <gdal/gdal.h>
#include <gdal/cpl_conv.h>
int main()
{
//zmienne
GDALDatasetH map_dataset;
GDALRasterBandH map_band;
GDALDriverH driver;
double geo_transform[6] ={0}; //lista informacji o geotransformacji (wspolrzedne pkt X = 0, Y = 0 na mapie), x_right punkt max na prawo, y_down punkt maks w dol
float *buff;
char *map_name, *points_name;
int i, cols, rows, nodata, pbSuccess; // i jest numeratorem,zmienna cols to ilosc kolumn pixeli mapy, rows ilosc wierszy pixeli mapy
float array_PCS[200];
int array_PIX[200];
GDALAllRegister();
map_name = (char*)malloc(101);
//wczytanie mapy
printf("Podaj nazwe pliku (mapa): ");
scanf("%100s", map_name);
map_dataset = GDALOpen(map_name, GA_ReadOnly);
if (map_dataset == NULL)
{
printf("Nie moge otworzyc pliku (mapa) \n");
free(map_name);
exit(0);
}
//wczytywanie sterownika
driver = GDALGetDriverByName("GTiff");
if(driver==NULL)
{
printf("blad sterownika\n");
GDALClose(map_dataset);
free(map_name);
exit(0);
}
map_band = GDALGetRasterBand(map_dataset, 1);
cols = GDALGetRasterBandXSize(map_band);
rows = GDALGetRasterBandYSize(map_band);
nodata = GDALGetRasterNoDataValue(map_band, &pbSuccess);
GDALGetGeoTransform(map_dataset, geo_transform);
//wczytanie pliku ze wspolrzednymi punktow
points_name = (char*)malloc(101);
printf("Podaj nazwe pliku (punkty): ");
scanf("%100s", points_name);
FILE* points_file = fopen(points_name, "r");
if (points_name == NULL)
{
printf("Nie moge otworzyc pliku (punkty)\n");
GDALClose(map_dataset);
free(map_name);
free(points_name);
exit(0);
}
while (fgetc(points_file)!='\n'); // pomija pierwsza linie pliku punkty.txt bo ta zawiera X Y a nie liczby
for (i=0; i<200; i+=2) //petla czytajaca wartosci z pliku i wpusujaca je do atblicy array_XY
{
if (feof(points_file)) //sprawdza czy petrla dotarla do konca pliku
{
break; //jezeli dotarla to wychodzi z petli
}
fscanf(points_file, "%f %f", &array_PCS[i], &array_PCS[i+1]); //tutaj zapisuje wpolrzedne plasie na liste array_PCS
array_PIX[i] = (int)(array_PCS[i] - geo_transform[0] - array_PCS[i+1] * geo_transform[2]) / geo_transform[1]; // konwertuje wspolrzedne plaskie X na wartosc w pixelach i wpisuje na liste array_PIX
array_PIX[i+1] = (int)(array_PCS[i+1] - geo_transform[3] - array_PCS[i] * geo_transform[4]) / geo_transform[5];
}
FILE* wynik = fopen("wynik.txt", "w"); //tworzy plik wyjsciowy o nazwie wynik.txt
buff = (float *)CPLMalloc(1* sizeof(float));
for (int j=0; j<i-3; j+=2) //petla w ktorej szuka pkt na warstwie i pobiera ich wysokosc, w tym samym czasie zapisuje do pliku wynnik txt
{
if (array_PIX[j] >= 0 && array_PIX[j] <= cols && array_PIX[j+1] >= 0 && array_PIX[j+1] <= rows) //warunek sprawdzajacy czy interesujacy nas pkt jest w granicach warstwy
{
CPLErr err;
CPLErr err2;
err = GDALRasterIO(map_band, GF_Read, array_PIX[j], array_PIX[j+1], 1, 1, buff, 1, 1, GDT_Float32, 0, 0);
err2 = GDALRasterIO(map_band, GF_Read, array_PIX[j+2], array_PIX[j+3], 1, 1, buff, 1, 1, GDT_Float32, 0, 0);
if(err>CE_Warning)
{
printf("Blad bufforu\n"); //wyswietla tekst przy bledzie
break; //jesli blad to wychodzi z petli
}
fprintf(wynik, "X: %.3f Y: %.3f", array_PCS[j], array_PCS[j+1]); //przepisuje wspolrzedne z listy
fprintf(wynik, "X: %.3f Y: %.3f", array_PCS[j+2], array_PCS[j+3]);
if (buff[0] == nodata) // sprawdza czy punkt ma wartosc wysokosci
{
fprintf(wynik, "Brak danych wysokosciowych");
}
else
{
fprintf(wynik, "Wysokosc: %.2f", buff[0]); //to dokleja wartosc wysokosci
}
fprintf(wynik, "\n");
}
else
{
fprintf(wynik, "X: %.3f Y: %.3f Punkt poza mapa\n", array_PCS[j], array_PCS[j+1]);
}
}
GDALClose(map_dataset);
free(map_name);
free(points_name);
CPLFree(buff);
}
Tutaj opis jak powinien działać ten program:
Program odczytuje wielkości z mapy rastrowej wzdłuż zadanego profilu. Program wczytuje
plik w formacie obsługiwanym przez GDAL (np. GeoTIFF). Nazwy plików wejściowego i
wynikowego oraz współrzędne podawane są jako parametry uruchomienia programu lub, w
przypadku nie podania parametrów, wczytywane interaktywnie. Podawane są współrzędne x,
y punktów rozpoczynającego i kończącego profil. Profil jest odcinkiem prostym. Wynikiem
działania programu jest lista współrzędnych x, y komórek wzdłuż profilu oraz wartości z
odczytane z warstwy wejściowej. Wynik zapisywany jest do pliku tekstowego. Przykładowymi
danymi jest fragmentu modelu DTM.
Będę wdzięczny za wszelkie sugestie.