Skip to content

Commit e1907bb

Browse files
rijnbrijnb
committed
Added boundaries check
1 parent d419cac commit e1907bb

File tree

2 files changed

+128
-30
lines changed

2 files changed

+128
-30
lines changed

‎mapcode‎

-696 KB
Binary file not shown.

‎mapcode.c‎

Lines changed: 128 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99

1010
static const double PI = 3.14159265358979323846;
1111
static const int RESULTS_MAX = 64;
12+
static const int SHOW_PROGRESS = 250;
1213

1314
static void usage(const char* appName) {
1415
printf("MAPCODE (C library version %s)\n", mapcode_cversion);
@@ -25,11 +26,13 @@ static void usage(const char* appName) {
2526
printf(" Encode a lat/lon to a Mapcode. If the territory code is specified, the\n");
2627
printf(" encoding will only succeeed if the lat/lon is located in the territory.\n");
2728
printf("\n");
28-
printf(" %s [-g | --grid] <nrPoints> [<seed>]\n", appName);
29+
printf(" %s [-b | --boundaries]\n", appName);
30+
printf(" %s [-g | --grid] <nrPoints>\n", appName);
2931
printf(" %s [-r | --random] <nrPoints> [<seed>]\n", appName);
3032
printf("\n");
31-
printf(" Create a test set of a number of a grid or random uniformly distributed\n");
32-
printf(" lat/lon pairs, (x, y, z) points and the Mapcode aliases.\n");
33+
printf(" Create a test set of lat/lon pairs based on the Mapcode boundaries database\n");
34+
printf(" as a fixed 3D grid or random uniformly distributed set of lat/lons with their\n");
35+
printf(" (x, y, z) coordinates and all Mapcode aliases.\n");
3336
printf("\n");
3437
printf(" The output format is:\n");
3538
printf(" <number-of-aliases> <lat-deg> <lon-deg> <x> <y> <z>\n");
@@ -46,30 +49,67 @@ static void usage(const char* appName) {
4649
printf(" and the (x, y, z) coordinates are placed on a sphere with radius 1.\n");
4750
}
4851

52+
static double radToDeg(double rad){
53+
return (rad / PI) * 180.0;
54+
}
55+
56+
static double degToRad(double deg){
57+
return (deg / 180.0) * PI;
58+
}
59+
4960
/**
5061
* Given a single number between 0..1, generate a latitude, longitude (in degrees) and a 3D
5162
* (x, y, z) point on a sphere with a radius of 1.
5263
*/
53-
static void unitToLatLonDegXYZ(
54-
const double unit1, const double unit2,
55-
double* latDeg, double* lonDeg, double* x, double* y, double* z) {
64+
static void unitToLatLonDeg(
65+
const double unit1, const double unit2, double* latDeg, double* lonDeg) {
5666

5767
// Calculate uniformly distributed 3D point on sphere (radius = 1.0):
5868
// http://mathproofs.blogspot.co.il/2005/04/uniform-random-distribution-on-sphere.html
5969
const double theta0 = (2.0 * PI) * unit1;
6070
const double theta1 = acos(1.0 - (2.0 * unit2));
61-
*x = sin(theta0) * sin(theta1);
62-
*y = cos(theta0) * sin(theta1);
63-
*z = cos(theta1);
71+
double x = sin(theta0) * sin(theta1);
72+
double y = cos(theta0) * sin(theta1);
73+
double z = cos(theta1);
6474

6575
// Convert Carthesian 3D point into lat/lon (radius = 1.0):
6676
// http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
67-
const double latRad = asin(*z);
68-
const double lonRad = atan2(*y, *x);
77+
const double latRad = asin(z);
78+
const double lonRad = atan2(y, x);
6979

7080
// Convert radians to degrees.
71-
*latDeg = isnan(latRad) ? 90.0 : (latRad * (180.0 / PI));
72-
*lonDeg = isnan(lonRad) ? 180.0 : (lonRad * (180.0 / PI));
81+
*latDeg = isnan(latRad) ? 90.0 : radToDeg(latRad);
82+
*lonDeg = isnan(lonRad) ? 180.0 : radToDeg(lonRad);
83+
}
84+
85+
static void convertLatLonToXYZ(double latDeg, double lonDeg, double* x, double* y, double* z) {
86+
double latRad = degToRad(latDeg);
87+
double lonRad = degToRad(lonDeg);
88+
*x = cos(latRad) * cos(lonRad);
89+
*y = cos(latRad) * sin(lonRad);
90+
*z = sin(latRad);
91+
}
92+
93+
static int printMapcodes(double lat, double lon, int iShowError) {
94+
const char* results[RESULTS_MAX];
95+
int context = 0;
96+
const int nrResults = coord2mc(results, lat, lon, context);
97+
if (nrResults <= 0) {
98+
if (iShowError) {
99+
fprintf(stderr, "error: cannot encode lat=%f, lon=%f)\n", lat, lon);
100+
}
101+
return -1;
102+
}
103+
double x;
104+
double y;
105+
double z;
106+
convertLatLonToXYZ(lat, lon, &x, &y, &z);
107+
printf("%d %lf %lf %lf %lf %lf\n", nrResults, lat, lon, x, y, z);
108+
for (int j = 0; j < nrResults; ++j) {
109+
printf("%s %s\n", results[(j * 2) + 1], results[(j * 2)]);
110+
}
111+
printf("\n");
112+
return nrResults;
73113
}
74114

75115
int main(const int argc, const char** argv)
@@ -135,6 +175,76 @@ int main(const int argc, const char** argv)
135175
printf("%s %s\n", results[(i * 2) + 1], results[(i * 2)]);
136176
}
137177
}
178+
else if ((strcmp(cmd, "-b") == 0) || (strcmp(cmd, "--boundaries") == 0)) {
179+
180+
// ------------------------------------------------------------------
181+
// Generate a test set based on the Mapcode boundaries.
182+
// ------------------------------------------------------------------
183+
if (argc != 2) {
184+
fprintf(stderr, "error: incorrect number of arguments\n\n");
185+
usage(appName);
186+
return -1;
187+
}
188+
for (int i = 0; i < NR_RECS; ++i) {
189+
long minLonE6;
190+
long maxLonE6;
191+
long minLatE6;
192+
long maxLatE6;
193+
double minLon;
194+
double maxLon;
195+
double minLat;
196+
double maxLat;
197+
double lat;
198+
double lon;
199+
200+
get_boundaries(i, &minLonE6, &minLatE6, &maxLonE6, &maxLatE6);
201+
minLon = ((double) minLonE6) / 1.0E6;
202+
maxLon = ((double) maxLonE6) / 1.0E6;
203+
minLat = ((double) minLatE6) / 1.0E6;
204+
maxLat = ((double) maxLatE6) / 1.0E6;
205+
206+
// Try center.
207+
lat = (maxLat - minLat ) / 2.0;
208+
lon = (maxLon - minLon ) / 2.0;
209+
210+
// Try center.
211+
printMapcodes(lat, lon, 0);
212+
213+
// Try corners.
214+
printMapcodes(minLat, minLon, 0);
215+
printMapcodes(minLat, maxLon, 0);
216+
printMapcodes(maxLat, minLon, 0);
217+
printMapcodes(maxLat, maxLon, 0);
218+
219+
// Try JUST inside.
220+
double factor = 1.0;
221+
for (int j = 1; j < 6; ++j) {
222+
223+
double d = 1.0 / factor;
224+
printMapcodes(minLat + d, minLon + d, 0);
225+
printMapcodes(minLat + d, maxLon - d, 0);
226+
printMapcodes(maxLat - d, minLon + d, 0);
227+
printMapcodes(maxLat - d, maxLon - d, 0);
228+
229+
// Try JUST outside.
230+
printMapcodes(minLat - d, minLon - d, 0);
231+
printMapcodes(minLat - d, maxLon + d, 0);
232+
printMapcodes(maxLat + d, minLon - d, 0);
233+
printMapcodes(maxLat + d, maxLon + d, 0);
234+
factor = factor * 10.0;
235+
}
236+
237+
// Try 22m outside.
238+
printMapcodes(minLat - 22, (maxLon - minLon) / 2, 0);
239+
printMapcodes(minLat - 22, (maxLon - minLon) / 2, 0);
240+
printMapcodes(maxLat + 22, (maxLon - minLon) / 2, 0);
241+
printMapcodes(maxLat + 22, (maxLon - minLon) / 2, 0);
242+
243+
if ((i % SHOW_PROGRESS) == 0) {
244+
fprintf(stderr, "Processed %d of %d regions...\r", i, NR_RECS);
245+
}
246+
}
247+
}
138248
else if ((strcmp(cmd, "-g") == 0) || (strcmp(cmd, "--grid") == 0) ||
139249
(strcmp(cmd, "-r") == 0) || (strcmp(cmd, "--random") == 0)) {
140250

@@ -177,18 +287,12 @@ int main(const int argc, const char** argv)
177287
double lonLargestNrOfResults = 0.0;
178288
int totalNrOfResults = 0;
179289

180-
const char* results[RESULTS_MAX];
181-
int context = 0;
182-
183290
int gridX = 0;
184291
int gridY = 0;
185292
int line = round(sqrt(nrPoints));
186293
for (int i = 0; i < nrPoints; ++i) {
187294
double lat;
188295
double lon;
189-
double x;
190-
double y;
191-
double z;
192296
double unit1;
193297
double unit2;
194298

@@ -209,25 +313,19 @@ int main(const int argc, const char** argv)
209313
}
210314
}
211315

212-
unitToLatLonDegXYZ(unit1, unit2, &lat, &lon, &x, &y, &z);
213-
const int nrResults = coord2mc(results, lat, lon, context);
214-
if (nrResults <= 0) {
215-
fprintf(stderr, "error: cannot encode lat=%f, lon=%f)\n", lat, lon);
216-
return -1;
217-
}
316+
unitToLatLonDeg(unit1, unit2, &lat, &lon);
317+
const int nrResults = printMapcodes(lat, lon, 1);
218318
if (nrResults > largestNrOfResults) {
219319
largestNrOfResults = nrResults;
220320
latLargestNrOfResults = lat;
221321
lonLargestNrOfResults = lon;
222322
}
223323
totalNrOfResults += nrResults;
224-
printf("%d %lf %lf %lf %lf %lf\n", nrResults, lat, lon, x, y, z);
225-
for (int j = 0; j < nrResults; ++j) {
226-
printf("%s %s\n", results[(j * 2) + 1], results[(j * 2)]);
324+
if ((i % SHOW_PROGRESS) == 0) {
325+
fprintf(stderr, "Created %d of %d 3D %s data points...\r", i, nrPoints, random ? "random" : "grid");
227326
}
228-
printf("\n");
229327
}
230-
fprintf(stderr, "Statistics:\n");
328+
fprintf(stderr, "\nStatistics:\n");
231329
fprintf(stderr, "Total number of 3D points generated = %d\n", nrPoints);
232330
fprintf(stderr, "Total number of Mapcodes generated = %d\n", totalNrOfResults);
233331
fprintf(stderr, "Average number of Mapcodes per 3D point = %f\n", ((float) totalNrOfResults) / ((float) nrPoints));

0 commit comments

Comments
 (0)