GRASS Programmer's Manual
6.4.2(2012)
|
00001 00017 #include <stdio.h> 00018 #include <string.h> 00019 #include <ctype.h> 00020 00021 #include <grass/gis.h> 00022 #include <grass/gprojects.h> 00023 #include <grass/glocale.h> 00024 00025 /* a couple defines to simplify reading the function */ 00026 #define MULTIPLY_LOOP(x,y,c,m) \ 00027 do {\ 00028 int i; \ 00029 for (i = 0; i < c; ++i) {\ 00030 x[i] *= m; \ 00031 y[i] *= m; \ 00032 }\ 00033 } while (0) 00034 00035 #define DIVIDE_LOOP(x,y,c,m) \ 00036 do {\ 00037 int i; \ 00038 for (i = 0; i < c; ++i) {\ 00039 x[i] /= m; \ 00040 y[i] /= m; \ 00041 }\ 00042 } while (0) 00043 00044 static double METERS_in = 1.0, METERS_out = 1.0; 00045 00063 int pj_do_proj(double *x, double *y, 00064 struct pj_info *info_in, struct pj_info *info_out) 00065 { 00066 int ok; 00067 double u, v; 00068 double h = 0.0; 00069 00070 METERS_in = info_in->meters; 00071 METERS_out = info_out->meters; 00072 00073 if (strncmp(info_in->proj, "ll", 2) == 0) { 00074 if (strncmp(info_out->proj, "ll", 2) == 0) { 00075 u = (*x) / RAD_TO_DEG; 00076 v = (*y) / RAD_TO_DEG; 00077 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h); 00078 *x = u * RAD_TO_DEG; 00079 *y = v * RAD_TO_DEG; 00080 } 00081 else { 00082 u = (*x) / RAD_TO_DEG; 00083 v = (*y) / RAD_TO_DEG; 00084 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h); 00085 *x = u / METERS_out; 00086 *y = v / METERS_out; 00087 } 00088 } 00089 else { 00090 if (strncmp(info_out->proj, "ll", 2) == 0) { 00091 u = *x * METERS_in; 00092 v = *y * METERS_in; 00093 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h); 00094 *x = u * RAD_TO_DEG; 00095 *y = v * RAD_TO_DEG; 00096 } 00097 else { 00098 u = *x * METERS_in; 00099 v = *y * METERS_in; 00100 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h); 00101 *x = u / METERS_out; 00102 *y = v / METERS_out; 00103 } 00104 } 00105 if (ok < 0) { 00106 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok)); 00107 } 00108 return ok; 00109 } 00110 00134 int pj_do_transform(int count, double *x, double *y, double *h, 00135 struct pj_info *info_in, struct pj_info *info_out) 00136 { 00137 int ok; 00138 int has_h = 1; 00139 00140 METERS_in = info_in->meters; 00141 METERS_out = info_out->meters; 00142 00143 if (h == NULL) { 00144 int i; 00145 00146 h = G_malloc(sizeof *h * count); 00147 /* they say memset is only guaranteed for chars ;-( */ 00148 for (i = 0; i < count; ++i) 00149 h[i] = 0.0; 00150 has_h = 0; 00151 } 00152 if (strncmp(info_in->proj, "ll", 2) == 0) { 00153 if (strncmp(info_out->proj, "ll", 2) == 0) { 00154 DIVIDE_LOOP(x, y, count, RAD_TO_DEG); 00155 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h); 00156 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG); 00157 } 00158 else { 00159 DIVIDE_LOOP(x, y, count, RAD_TO_DEG); 00160 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h); 00161 DIVIDE_LOOP(x, y, count, METERS_out); 00162 } 00163 } 00164 else { 00165 if (strncmp(info_out->proj, "ll", 2) == 0) { 00166 MULTIPLY_LOOP(x, y, count, METERS_in); 00167 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h); 00168 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG); 00169 } 00170 else { 00171 MULTIPLY_LOOP(x, y, count, METERS_in); 00172 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h); 00173 DIVIDE_LOOP(x, y, count, METERS_out); 00174 } 00175 } 00176 if (!has_h) 00177 G_free(h); 00178 00179 if (ok < 0) { 00180 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok)); 00181 } 00182 return ok; 00183 }