Weter, on 05 September 2024 - 11:34 AM, said:
Thank you very much!
Posted 11 September 2024 - 05:23 PM
Posted 11 September 2024 - 06:12 PM
Posted 11 September 2024 - 06:38 PM
Posted 11 September 2024 - 06:48 PM
Posted 11 September 2024 - 07:35 PM
Posted 12 September 2024 - 04:58 AM
#include <gdal/gdal.h> float getElevation(double lat, double lng, float* elevData, int xSize, int ySize, double* transform) { double x= (lng-transform[0])/transform[1]-.5; double y= (lat-transform[3])/transform[5]-.5; if (x<0 || y<0 || x>=xSize-1 || y>=ySize-1) { // fprintf(stderr,"out of range %f %f %f %f\n",x,y,lat,lng); return -1; } int x0= (int)floor(x); int y0= (int)floor(y); int x1= x0+1; int y1= y0+1; double xa= x-x0; double ya= y-y0; float a00= elevData[x0+y0*xSize]; float a01= elevData[x1+y0*xSize]; float a10= elevData[x0+y1*xSize]; float a11= elevData[x1+y1*xSize]; return (1-ya)*((1-xa)*a00+xa*a01) + ya*((1-xa)*a10+xa*a11); } void saveTerrain() { GDALAllRegister(); GDALDatasetH rasterData= GDALOpen(rasterFile.c_str(),GA_ReadOnly); if (rasterData == NULL) throw "cannot read raster file"; int xSize= GDALGetRasterXSize(rasterData); int ySize= GDALGetRasterYSize(rasterData); fprintf(stderr,"size %d %d %d\n",xSize,ySize, GDALGetRasterCount(rasterData)); double transform[6]; GDALGetGeoTransform(rasterData,transform); fprintf(stderr,"transform"); for (int i=0; i<6; i++) fprintf(stderr," %f",transform[i]); fprintf(stderr,"\n"); GDALRasterBandH band= GDALGetRasterBand(rasterData,1); fprintf(stderr,"band size %d\n",GDALGetRasterBandXSize(band)); float* elevData= (float*) malloc(xSize*ySize*sizeof(float)); if (elevData == NULL) throw "cannot allocate memory"; if (GDALRasterIO(band,GF_Read,0,0,xSize,ySize,elevData,xSize,ySize, GDT_Float32,0,0)) throw "cannot read elevData\n"; for (MSTSRoute::TileMap::iterator i=mstsRoute->tileMap.begin(); i!=mstsRoute->tileMap.end(); ++i) { MSTSRoute::Tile* t= i->second; fprintf(stderr,"tile %d %d %f %f %f\n", t->x,t->z,t->floor,t->scale,t->floor+65535*t->scale); mstsRoute->readTerrain(t); float x0= 2048*(t->x-mstsRoute->centerTX); float z0= 2048*(t->z-mstsRoute->centerTZ); fprintf(stderr," %f %f\n",x0,z0); float minA= 1e10; float maxA= 0; for (int j=0; j<256; j++) { for (int k=0; k<256; k++) { double x= x0+8*(k-128); double z= z0+8*(128-j); double lng,lat; xy2ll(x,z,&lat,&lng); float a= getElevation(lat,lng, elevData,xSize,ySize,transform); #if 0 if (j%128==0 && k%128==0) fprintf(stderr, "ll %d %d %d %d %f %f %f %f %f\n", t->x,t->z,j,k,lat,lng,a, (lng-transform[0])/transform[1]-.5, (lat-transform[3])/transform[5]-.5); #endif if (a < 0) a= 0; if (minA > a) minA= a; if (maxA < a) maxA= a; } } fprintf(stderr,"minmax %f %f\n",minA,maxA); saveFloorScale(t,minA-1,(maxA-minA+2)/65535); for (int j=0; j<256; j++) { for (int k=0; k<256; k++) { double x= x0+8*(k-128); double z= z0+8*(128-j); double lng,lat; xy2ll(x,z,&lat,&lng); float a= getElevation(lat,lng, elevData,xSize,ySize,transform); if (a < 0) a= 0; t->terrain->y[j][k]= a<t->floor ? 0 : a>65535*t->scale+t->floor ? 65535 : (int)((a-t->floor)/t->scale); // fprintf(stderr,"%d %d %d %f\n", // j,k,t->terrain->y[j][k],a); } } mstsRoute->writeTerrain(t); t->freeTerrain(); } free(elevData); GDALClose(rasterData); }
Posted 12 September 2024 - 03:32 PM
eric from trainsim, on 11 September 2024 - 07:35 PM, said: