pd.if.org Git - isea/blob - pjexamples.c

   1 /*
   2  * options
   3  *
   4  * orientation
   5  * -l longitude
   6  * -z azimuth
   7  * -p latitude
   8  *
   9  * -I isea standard orientation
  10  * -D dymaxion orientation
  11  * -P vertex at pole
  12  *
  13  * -a aperture
  14  * -r resolution
  15  * -t h t d hex tri diamond
  16  *
  17  *  addressing
  18  *  -O q2dd seqnum interleave plane q2dd projtri vertex2dd
  19  *  -I
  20  *
  21  * operation
  22  * -g generate grid
  23  * -x transform
  24  * -b bin point values (average)
  25  * -e bin point presence
  26  * -G spacing resolution generate graticule
  27  */
  28
  29 #if 0
  30 void usage(void) {
  31         fprintf(stderr, "%s", "isea [-l <lon>] [-z <azimuth>] [-p <latitude>] [-I] [-D] [-P] [-a <aperture>] [-r <resolution>] [-t <topology>] [-A <addressing>] [-g] [-b] [-e] [-x]\n");
  32 }
  33
  34 void die(char *s) {
  35         fprintf(stderr, "%s\n", s);
  36         exit(EXIT_FAILURE);
  37 }
  38
  39 static int verbose = 0;
  40 static int precision = 10;
  41
  42 void options(int ac, char **av, struct isea_dgg *g) {
  43         int i = 1;
  44         char *optarg;
  45
  46         for (i=1;i<ac;i++) {
  47                 if (av[i][0] != '-') {
  48                         break;
  49                 }
  50                 switch (av[i][1]) {
  51                         case 'v':
  52                                 verbose = 1; break;
  53                         case 'I':
  54                                 isea_orient_isea(g); break;
  55                         case 'P':
  56                                 isea_orient_pole(g); break;
  57                         case 'D':
  58                                 isea_orient_dymax(g); break;
  59                         case 'R':
  60                                 g->radius = strtod(av[++i],NULL); break;
  61                         case 'l':
  62                                 g->o_lon = strtod(av[++i],NULL); break;
  63                         case 'p':
  64                                 g->o_lat = strtod(av[++i],NULL); break;
  65                         case 'z':
  66                                 g->o_az = strtod(av[++i],NULL); break;
  67                         case 'd':
  68                                 precision = atoi(av[++i]); break;
  69                         case 'a':
  70                                 g->aperture = atoi(av[++i]); break;
  71                         case 'r':
  72                                 g->resolution = atoi(av[++i]); break;
  73                         case 't':
  74                                 switch (av[++i][0]) {
  75                                         case 'h':
  76                                                 g->topology = ISEA_HEXAGON;
  77                                                 break;
  78                                         case 't':
  79                                                 g->topology = ISEA_TRIANGLE;
  80                                                 break;
  81                                         case 'd':
  82                                                 g->topology = ISEA_DIAMOND;
  83                                                 break;
  84                                         default:
  85                                                 usage();
  86                                                 die("Invalid topology");
  87                                                 break;
  88                                 };
  89                                 break;
  90                         case 'O':
  91                                 if (av[i][2]) {
  92                                         optarg = &av[i][2];
  93                                 } else {
  94                                         optarg = av[++i];
  95                                 }
  96                                 if (!strcmp(optarg, "projtri")) {
  97                                         g->output = ISEA_PROJTRI;
  98                                 }
  99                                 else if (!strcmp(optarg, "vertex2dd")) {
 100                                         g->output = ISEA_VERTEX2DD;
 101                                 }
 102                                 else if (!strcmp(optarg, "q2dd")) {
 103                                         g->output = ISEA_Q2DD;
 104                                 }
 105                                 else if (!strcmp(optarg, "q2di")) {
 106                                         g->output = ISEA_Q2DI;
 107                                 }
 108                                 else if (!strcmp(optarg, "seqnum")) {
 109                                         g->output = ISEA_SEQNUM;
 110                                 }
 111                                 else if (!strcmp(optarg, "plane")) {
 112                                         g->output = ISEA_PLANE;
 113                                 } else {
 114                                         fprintf(stderr, "unknown output format: %s\n",optarg);
 115                                         die("exiting");
 116                                 }
 117                                 break;
 118                         default:
 119                                 usage();
 120                                 die("Unknown option");
 121                 }
 122         }
 123 }
 124
 125 #endif
 126
 127 #define PROJ_PARMS__ \
 128         struct isea_dgg dgg;
 129
 130 #define PJ_LIB__
 131 #include <projects.h>
 132
 133 PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
 134
 135 FORWARD(s_forward);
 136         struct isea_pt out;
 137         struct isea_geo in;
 138
 139         in.lon = lp.lam;
 140         in.lat = lp.phi;
 141
 142         out = isea_forward(&P->dgg, &in);
 143
 144         xy.x = out.x;
 145         xy.y = out.y;
 146
 147         return xy;
 148 }
 149 FREEUP; if (P) pj_dalloc(P); }
 150
 151 ENTRY0(isea)
 152         char *opt;
 153
 154         P->fwd = s_forward;
 155         isea_grid_init(&P->dgg);
 156
 157         P->dgg.output = ISEA_PLANE;
 158 /*              P->dgg.radius = P->a; /* otherwise defaults to 1 */
 159         /* calling library will scale, I think */
 160
 161         opt = pj_param(P->params, "sorient").s;
 162         if (opt) {
 163                 if (!strcmp(opt, "isea")) {
 164                         isea_orient_isea(&P->dgg);
 165                 } else if (!strcmp(opt, "pole")) {
 166                         isea_orient_pole(&P->dgg);
 167                 } else {
 168                         E_ERROR(-34);
 169                 }
 170         }
 171
 172         if (pj_param(P->params, "tazi").i) {
 173                 P->dgg.o_az = pj_param(P->params, "razi").f;
 174         }
 175
 176         if (pj_param(P->params, "tlon_0").i) {
 177                 P->dgg.o_lon = pj_param(P->params, "rlon_0").f;
 178         }
 179
 180         if (pj_param(P->params, "tlat_0").i) {
 181                 P->dgg.o_lat = pj_param(P->params, "rlat_0").f;
 182         }
 183
 184         if (pj_param(P->params, "taperture").i) {
 185                 P->dgg.aperture = pj_param(P->params, "iaperture").i;
 186         }
 187
 188         if (pj_param(P->params, "tresolution").i) {
 189                 P->dgg.resolution = pj_param(P->params, "iresolution").i;
 190         }
 191
 192         opt = pj_param(P->params, "smode").s;
 193         if (opt) {
 194                 if (!strcmp(opt, "plane")) {
 195                         P->dgg.output = ISEA_PLANE;
 196                 } else if (!strcmp(opt, "di")) {
 197                         P->dgg.output = ISEA_Q2DI;
 198                 }
 199                 else if (!strcmp(opt, "dd")) {
 200                         P->dgg.output = ISEA_Q2DD;
 201                 }
 202                 else if (!strcmp(opt, "hex")) {
 203                         P->dgg.output = ISEA_HEX;
 204                 }
 205                 else {
 206                         E_ERROR(-34);
 207                 }
 208         }
 209
 210         if (pj_param(P->params, "trescale").i) {
 211                 P->dgg.radius = ISEA_SCALE;
 212         }
 213
 214         if (pj_param(P->params, "tresolution").i) {
 215                 P->dgg.resolution = pj_param(P->params, "iresolution").i;
 216         } else {
 217                 P->dgg.resolution = 4;
 218         }
 219
 220         if (pj_param(P->params, "taperture").i) {
 221                 P->dgg.aperture = pj_param(P->params, "iaperture").i;
 222         } else {
 223                 P->dgg.aperture = 3;
 224         }
 225
 226 ENDENTRY(P)
 227
 228 #if 0
 229 int
 230 main(int ac, char *av[])
 231 {
 232         struct isea_dgg dgg;
 233         struct isea_geo in;
 234         struct isea_pt  out, coord;
 235         int             tri,calctri;
 236         int             graticule = 0;
 237         int     quad = 0;
 238
 239         char            line[1024];
 240         char           *rest, *nl;
 241
 242         isea_grid_init(&dgg);
 243         dgg.output = ISEA_PLANE;
 244         dgg.radius = ISEA_SCALE; /* 0.8301572857837594396028083; */
 245
 246         options(ac, av, &dgg);
 247
 248         if (0) {
 249                 printf("isea lat = %lf, lon = %lf, az = %lf\n",
 250                        dgg.o_lat * 180.0 / M_PI, dgg.o_lon * 180.0 / M_PI,
 251                        dgg.o_az * 180.0 / M_PI);
 252         }
 253
 254         while (fgets(line, sizeof line, stdin)) {
 255                 line[sizeof line-1] = 0;
 256                 if (line[0] == '#') {
 257                         continue;
 258                 }
 259                 if (!strncmp(line, "ENDPOLY", 7)) {
 260                         printf("%s", line);
 261                         continue;
 262                 }
 263                 in.lon = strtod(line, &rest);
 264                 in.lat = strtod(rest, &rest);
 265
 266                 nl = strchr(rest, '\n');
 267                 if (nl) *nl = 0;
 268
 269                 in.lon *= M_PI / 180.0;
 270                 in.lat *= M_PI / 180.0;
 271
 272                 out = isea_forward(&dgg, &in);
 273                 tri = dgg.triangle;
 274                 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
 275                 coord = out;
 276
 277                 if (dgg.output == ISEA_PLANE) {
 278                         calctri = isea_xy_triangle(&coord,
 279                                 dgg.output == ISEA_PLANE ? dgg.radius : 1.0);
 280                         if (tri != calctri) {
 281                                 fprintf(stderr,
 282                                         "got %d (%f, %f), expected %d\n",
 283                                         calctri, out.x, out.y, tri);
 284                         }
 285                 }
 286
 287                 in.lon *= 180.0 / M_PI;
 288                 in.lat *= 180.0 / M_PI;
 289
 290                 switch (dgg.output) {
 291                 case ISEA_PLANE:
 292                         printf("%d %.*f %.*f%s\n", tri, precision, out.x, precision, out.y, rest);
 293                         break;
 294                 case ISEA_PROJTRI:
 295                         printf("%d %.*f %.*f%s\n", tri-1, precision, out.x,
 296                                precision, out.y, rest);
 297                         break;
 298                 case ISEA_VERTEX2DD:
 299                         printf("%d %d %.*f %.*f%s\n", quad, tri-1, precision,
 300                                         out.x, precision, out.y, rest);
 301                         break;
 302                 case ISEA_Q2DD:
 303                         /* Same as above, we just don't print as much */
 304                         printf("%d %.*f %.*f%s\n", quad, precision,
 305                                         out.x, precision, out.y, rest);
 306                         break;
 307                 case ISEA_Q2DI:
 308                         printf("%d %.0lf %.0lf%s\n", dgg.quad, coord.x, coord.y, rest);
 309                         break;
 310                 case ISEA_SEQNUM:
 311                         quad = isea_ptdi(&dgg, tri, &out, &coord);
 312                         printf("%d%s\n", isea_disn(&dgg, quad, &coord), rest);
 313                         break;
 314                 }
 315                 fflush(stdout);
 316         }
 317
 318         if (graticule) {
 319
 320         }
 321
 322         exit(EXIT_SUCCESS);
 323 }
 324 #endif