From b7a8170d4ae7bf30463b3ade7f3c004d4d8fe062 Mon Sep 17 00:00:00 2001 From: alex Date: Thu, 29 Sep 2016 22:19:50 +0100 Subject: [PATCH] 2.92 Bigger, faster LIDAR --- CHANGELOG | 4 ++ inputs.cc | 204 +++++++++++++++++------------------------------------ main.cc | 20 +++--- outputs.cc | 10 +-- 4 files changed, 82 insertions(+), 156 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 5d203f1..4d14a02 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,5 +1,9 @@ SIGNAL SERVER CHANGELOG +2.92 - 29 September 2016 +Overhauled LIDAR functions to read in (more) native format tiles. +Replaced 'super tile' code with more scalable multi-tile code: currently limited to 16. + 2.9 - 01 September 2016 MODIS Landcover support in ASCII Grid format. Works for 1200/3600 SRTM only. diff --git a/inputs.cc b/inputs.cc index 8307673..fc9660e 100644 --- a/inputs.cc +++ b/inputs.cc @@ -6,6 +6,7 @@ #include "common.h" #include "main.hh" + int loadClutter(char *filename, double radius, struct site tx) { /* This function reads a MODIS 17-class clutter file in ASCII Grid format. @@ -124,37 +125,29 @@ int loadClutter(char *filename, double radius, struct site tx) return 0; } -void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx, - double n, double e, double s, double west) +void readLIDAR(FILE *fd, int h, int w, int indx,double n, double e, double s, double west) { - int x = 0, y = 0, reads = 0; + int x = 0, y = 0, reads = 0, a=0, b=0, avg=0, tWidth = 0, tHeight = 0; char line[150000]; char *pch; - // use offsets to determine middle lat/lon for 5 x 10k tiles - // TALL - if (voffset == 0 && h == 10000) { - s = (s+n)/2; - h = 5000; - } - if (voffset == 5000 && h == 10000) { - n = (s+n)/2; - } - // WIDE - if (hoffset == 0 && w == 10000) { - e = (e+west)/2; - w = 5000; - } - if (hoffset == 5000 && w == 10000) { - west = (e+west)/2; - w = 5000; - } - dem[indx].max_north=n; dem[indx].min_west=e; dem[indx].min_north=s; dem[indx].max_west=west; + if (max_north == -90) + max_north = dem[indx].max_north; + + else if (dem[indx].max_north > max_north) + max_north = dem[indx].max_north; + + if (min_north == 90) + min_north = dem[indx].min_north; + + else if (dem[indx].min_north < min_north) + min_north = dem[indx].min_north; + if (dem[indx].max_west > max_west) max_west = dem[indx].max_west; if (dem[indx].min_west < min_west) @@ -186,18 +179,12 @@ void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx, for (y = h-1; y > -1; y--) { x = w-1; + if (fgets(line, 150000, fd) != NULL) { - // do nothing until n rows have passed - if (y < voffset || voffset == 0) { - pch = strtok(line, " "); - - //dummy reads until we reach offset - // for 5000 offset, width must be 10e3 - for (n = 0; n < hoffset; n++) - pch = strtok(NULL, " "); + pch = strtok(line, " "); // split line into values while (pch != NULL && x > -1) { - if (atoi(pch) < -999) + if (atoi(pch) <= -9999) pch = "0"; dem[indx].data[y][x] = atoi(pch); dem[indx].signal[x][y] = 0; @@ -210,46 +197,32 @@ void readLIDAR(FILE *fd, int hoffset, int voffset, int h, int w, int indx, dem[indx].min_el = atoi(pch); min_elevation = dem[indx].min_el; } - x--; + + x--; pch = strtok(NULL, " "); }//while - }//voffset + + } else { fprintf(stdout, "LIDAR error @ x %d y %d indx %d\n", x, y, indx); }//if }//for + } int loadLIDAR(char *filenames) { - /* - * This function reads either 9 LIDAR tiles of n rows and n columns - * in ASCII grid format OR a single super tile composed of 2 or more - * tiles. The tile must have WGS84 bounds in the header in the order: - * WEST,SOUTH,EAST,NORTH - * ncols 5000 - * nrows 5000 - * xllcorner -2.291359 - * yllcorner 51.788295 - * xurcorner -2.146674 - * yurcorner 51.878474 - * cellsize 2 - * NODATA_value -9999 - * - * Tiles must be entered in the format - * - * -lid tile1.asc,tile2.asc,tile3.asc - */ char *filename; - char *files[4]; // 4 tiles + char *files[25]; // 25 tiles int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, dem_alloced = 0; - double xll, yll, xur, yur, cellsize; + double xll, yll, xur, yur, cellsize, avgCellsize; char found, free_page = 0, jline[20], lid_file[255], path_plus_name[255], *junk = NULL; char line[50000]; char * pch; + double TO_DEG = (180 / PI); FILE *fd; // test for multiple files @@ -265,10 +238,10 @@ int loadLIDAR(char *filenames) if (fd != NULL) { - if (fgets(line, 19, fd) != NULL) { + if (fgets(line, 255, fd) != NULL) { pch = strtok (line," "); pch = strtok (NULL, " "); - width = atoi(pch); + width = atoi(pch); // ncols if (debug) { fprintf(stdout, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, width); @@ -276,63 +249,47 @@ int loadLIDAR(char *filenames) } if (!dem_alloced) { - IPPD = width; + IPPD = width * 1.1; // +10% ARRAYSIZE = (MAXPAGES * IPPD) + 10; do_allocs(); dem_alloced = 1; } } - if (fgets(line, 19, fd) != NULL) - height = atoi(pch); - fgets(line, 24, fd); + if (fgets(line, 255, fd) != NULL) + height = atoi(pch); // nrows - if (fgets(line, 24, fd) != NULL) { - //xll=atof(pch); - sscanf(pch, "%lf", &xll); - } - fgets(line, 24, fd); - if (fgets(line, 24, fd) != NULL) { - //yll=atof(pch); - sscanf(pch, "%lf", &yll); + if (fgets(line, 255, fd) != NULL) { + sscanf(pch, "%lf", &xll); // xll } - fgets(line, 24, fd); - - if (fgets(line, 24, fd) != NULL) { - //xur=atof(pch); - sscanf(pch, "%lf", &xur); + if (fgets(line, 255, fd) != NULL) { + sscanf(pch, "%lf", &yll); // yll } - fgets(line, 24, fd); + if (fgets(line, 255, fd) != NULL) + sscanf(pch, "%lf", &cellsize); + avgCellsize=avgCellsize+cellsize; - if (fgets(line, 24, fd) != NULL) { - //yur=atof(pch); - sscanf(pch, "%lf", &yur); + if(cellsize>=0.5){ // 50cm LIDAR? + // compute xur and yur with inverse haversine if cellsize in *metres* + double roundDistance = (width*cellsize)/6371000; + yur = asin(sin(yll*DEG2RAD) * cos(roundDistance) + cos(yll * DEG2RAD) * sin(roundDistance) * cos(0)) * TO_DEG; + xur = ((xll*DEG2RAD) + atan2(sin(90*DEG2RAD) * sin(roundDistance) * cos(yll*DEG2RAD), cos(roundDistance) - sin(yll * DEG2RAD) * sin(yur*DEG2RAD))) * TO_DEG; + }else{ + // Degrees with GDAL option: -co "FORCE_CELLSIZE=YES" + xur = xll+(cellsize*width); + yur = yll+(cellsize*height); } - - fgets(line, 15, fd); - - if (fgets(line, 15, fd) != NULL) - cellsize = strtod(pch, NULL); - - // LIDAR 10m @ 10800 PPD - if (cellsize == 10.0) - MAXRAD = 30; - // LIDAR 2m @ 54000 PPD - if (cellsize == 2.0) - MAXRAD = 15; - // LIDAR 1m @ 108000 PPD! - if (cellsize == 1.0) - MAXRAD = 10; - - if (xur < eastoffset) + + if (xur > eastoffset) eastoffset = xur; - if (xll > westoffset) + if (xll < westoffset) westoffset = xll; if (debug) - fprintf(stdout, "PRE yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta); + fprintf(stdout,"%d, %d, %.7f, %.7f, %.0f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur); + // Greenwich straddling hack if (xll < 0 && xur > 0) { @@ -353,58 +310,16 @@ int loadLIDAR(char *filenames) if (debug) fprintf(stdout, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta); - if (yll < min_north) - min_north = yll; - if (yur > max_north) - max_north = yur; - fgets(line, 30, fd); // NODATA + fgets(line, 255, fd); // NODATA pos = ftell(fd); // tile 0 [x| ] if (debug) - fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, 0, height, width, indx, yur, xur, yll, xll); + fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", height, width, indx, yur, xur, yll, xll); - readLIDAR(fd, 0, 0, height, width, indx, yur, xur, yll, xll); - //rewind - fseek(fd, pos, SEEK_SET); + readLIDAR(fd, height, width, indx, yur, xur, yll, xll); - // tile 1 [ |x] - if (width == 2000 || - width == 10000 || - width == 10082) { - indx++; - if (debug) - fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, 0, height, width, indx, yur, xur, yll, xll); - - readLIDAR(fd, width / 2, 0, height, width, indx, yur, xur, yll, xll); - } - //rewind - fseek(fd, pos, SEEK_SET); - - // tile 2 [x | ] - if (height == 2000 || - height == 10000 || - height == 10082) { - indx++; - if (debug) - fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", 0, height / 2, height, width, indx, yur, xur, yll, xll); - - readLIDAR(fd, 0, height / 2, height, width, indx, yur, xur, yll, xll); - } - //rewind - fseek(fd, pos, SEEK_SET); - - // tile 3 [ |x] - if ((width == 2000 && height == 2000) || - (width == 10000 && height == 10000) || - (width == 10082 && height == 10082)) { - indx++; - if (debug) - fprintf(stdout, "readLIDAR(fd,%d,%d,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", width / 2, height / 2, height, width, indx, yur, xur, yll, xll); - - readLIDAR(fd, width / 2, height / 2, height, width, indx, yur, xur, yll, xll); - } fclose(fd); if (debug) fprintf(stdout, "LIDAR LOADED %d x %d\n", width, height); @@ -414,6 +329,15 @@ int loadLIDAR(char *filenames) indx++; } // filename(s) + // TESTING! + IPPD=width; + ippd=width; + avgCellsize=avgCellsize/fc; + height = (unsigned)((max_north-min_north) / avgCellsize); + width = (unsigned)((max_west-min_west) / avgCellsize); + + if (debug) + fprintf(stdout, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west); return 0; } diff --git a/main.cc b/main.cc index d24a049..cca6bff 100644 --- a/main.cc +++ b/main.cc @@ -1,4 +1,4 @@ -double version = 2.91; +double version = 2.92; /****************************************************************************\ * Signal Server: Radio propagation simulator by Alex Farrant QCVS, 2E0TDW * ****************************************************************************** @@ -45,7 +45,7 @@ char string[255], sdf_path[255], udt_file[255], opened = 0, gpsav = double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd, fzone_clearance = 0.6, forced_freq, clutter, lat, lon, txh, tercon, terdic, north, east, south, west, dBm, loss, field_strength, - min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=-180, eastoffset=180, delta=0, rxGain=0; + min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=180, eastoffset=-180, delta=0, rxGain=0; int ippd, mpi, max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold, @@ -1024,7 +1024,7 @@ int main(int argc, char *argv[]) unsigned char LRmap = 0, txsites = 0, topomap = 0, geo = 0, kml = 0, area_mode = 0, max_txsites, ngs = 0; - char mapfile[255], udt_file[255], ano_filename[255], lidar_tiles[512], clutter_file[255]; + char mapfile[255], udt_file[255], ano_filename[255], lidar_tiles[1024], clutter_file[255]; double altitude = 0.0, altitudeLR = 0.0, tx_range = 0.0, rx_range = 0.0, deg_range = 0.0, deg_limit = 0.0, deg_range_lon; @@ -1036,9 +1036,9 @@ int main(int argc, char *argv[]) } if (strstr(argv[0], "signalserverLIDAR")) { - MAXPAGES = 4; + MAXPAGES = 25; // 5x5 lidar = 1; - IPPD = 5000; + IPPD = 500; // will be overridden based upon file header... } strncpy(ss_name, "Signal Server\0", 14); @@ -1227,7 +1227,7 @@ int main(int argc, char *argv[]) z = x + 1; lidar=1; if (z <= y && argv[z][0] && argv[z][0] != '-') - strncpy(lidar_tiles, argv[z], 253); + strncpy(lidar_tiles, argv[z], 1022); } if (strcmp(argv[x], "-res") == 0) { @@ -1616,14 +1616,18 @@ int main(int argc, char *argv[]) "correct and try again.\n"); exit(EXIT_FAILURE); } - ippd = IPPD; + if(debug){ - fprintf(stdout,"%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west); + fprintf(stdout,"%.4f,%.4f,%.4f,%.4f,%d x %d\n",max_north,min_west,min_north,max_west,width,height); } ppd=rint(height / (max_north-min_north)); yppd=rint(width / (max_west-min_west)); + //ppd=(double)ippd; + //yppd=ppd; + + if(delta>0){ tx_site[0].lon+=delta; } diff --git a/outputs.cc b/outputs.cc index 187a0c6..4e4ed95 100644 --- a/outputs.cc +++ b/outputs.cc @@ -567,18 +567,12 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, x0 = (int)rint((ppd * (lat - - (double)dem[indx].min_north))); // +4549 fix + (double)dem[indx].min_north))); y0 = mpi - - (int)rint(ppd * // ONLY exception. All others are yppd + (int)rint(ppd * (LonDiff ((double)dem[indx].max_west,lon))); - // fix for multi-tile lidar - if(width==10000 && (indx==1 || indx==3)){ - if(y0 >= 3432){ //3510,3535 - y0=y0-3432; - } - } if (x0 >= 0 && x0 <= mpi && y0 >= 0 && y0 <= mpi)