1 /****************************************************************************
2 *
3 * MODULE: r.latlong
4 * AUTHOR(S): Yann Chemin - yann.chemin@gmail.com
5 * PURPOSE: Calculates the longitude of the pixels in the map.
6 *
7 * COPYRIGHT: (C) 2002-2008, 2012 by the GRASS Development Team
8 *
9 * This program is free software under the GNU General
10 * Public License (>=v2). Read the file COPYING that
11 * comes with GRASS for details.
12 *
13 *****************************************************************************/
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <grass/gis.h>
19 #include <grass/raster.h>
20 #include <grass/gprojects.h>
21 #include <grass/glocale.h>
22
main(int argc,char * argv[])23 int main(int argc, char *argv[])
24 {
25 struct Cell_head cellhd; /*region+header info */
26 int nrows, ncols;
27 int row, col;
28 int not_ll = 0; /*if proj is not lat/long, it will be 1. */
29
30 struct GModule *module;
31 struct Option *input1, *output1;
32 struct Flag *flag1;
33 struct History history; /*metadata */
34 struct pj_info iproj, oproj, tproj;
35 struct Key_Value *in_proj_info, *in_unit_info;
36
37 /************************************/
38 char *result1; /*output raster name */
39
40 int infd;
41 int outfd1;
42 char *in;
43 double xmin, ymin;
44 double xmax, ymax;
45 double stepx, stepy;
46 double latitude, longitude;
47 void *inrast;
48 DCELL * outrast1;
49 DCELL d;
50
51 /************************************/
52 G_gisinit(argv[0]);
53
54 module = G_define_module();
55 G_add_keyword(_("raster"));
56 G_add_keyword(_("latitude"));
57 G_add_keyword(_("longitude"));
58 G_add_keyword(_("projection"));
59 module->description = _("Creates a latitude/longitude raster map.");
60
61 /* Define the different options */
62 input1 = G_define_standard_option(G_OPT_R_INPUT);
63
64 output1 = G_define_standard_option(G_OPT_R_OUTPUT);
65 output1->description = _("Name for output latitude or longitude raster map");
66
67 flag1 = G_define_flag();
68 flag1->key = 'l';
69 flag1->description = _("Longitude output");
70
71 /********************/
72 if (G_parser(argc, argv))
73 exit(EXIT_FAILURE);
74
75 in = input1->answer;
76 result1 = output1->answer;
77
78 /***************************************************/
79 infd = Rast_open_old(in, "");
80 Rast_get_cellhd(in, "", &cellhd);
81 inrast = Rast_allocate_d_buf();
82
83 /***************************************************/
84 xmin = cellhd.west;
85 xmax = cellhd.east;
86 ymin = cellhd.south;
87 ymax = cellhd.north;
88 nrows = Rast_window_rows();
89 ncols = Rast_window_cols();
90 stepx = abs(xmax-xmin)/(double)ncols;
91 stepy = abs(ymax-ymin)/(double)nrows;
92
93 /*Stolen from r.sun */
94 /* Set up parameters for projection to lat/long if necessary */
95 if ((G_projection() != PROJECTION_LL)) {
96 not_ll = 1;
97
98 if ((in_proj_info = G_get_projinfo()) == NULL)
99 G_fatal_error(_("Unable to get projection info of current location"));
100 if ((in_unit_info = G_get_projunits()) == NULL)
101 G_fatal_error(_("Unable to get projection units of current location"));
102 if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
103 G_fatal_error(_("Unable to get projection key values of current location"));
104 G_free_key_value(in_proj_info);
105 G_free_key_value(in_unit_info);
106
107 oproj.pj = NULL;
108 tproj.def = NULL;
109
110 if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
111 G_fatal_error(_("Unable to initialize coordinate transformation"));
112 } /* End of stolen from r.sun */
113
114 outrast1 = Rast_allocate_d_buf();
115
116 outfd1 = Rast_open_new(result1, DCELL_TYPE);
117
118 for (row = 0; row < nrows; row++)
119 {
120 G_percent(row, nrows, 2);
121
122 Rast_get_d_row(infd, inrast, row);
123
124 for (col = 0; col < ncols; col++)
125 {
126 latitude = ymax - ((double)row * stepy);
127 longitude = xmin + ((double)col * stepx);
128 if (not_ll) {
129 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
130 &longitude, &latitude, NULL) < 0)
131 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
132 "GPJ_transform()");
133 }
134 if (flag1->answer)
135 d = longitude;
136 else
137 d = latitude;
138 outrast1[col] = d;
139 }
140 Rast_put_d_row(outfd1, outrast1);
141 }
142 G_free(inrast);
143 Rast_close(infd);
144 G_free(outrast1);
145 Rast_close(outfd1);
146
147 Rast_short_history(result1, "raster", &history);
148 Rast_command_history(&history);
149 Rast_write_history(result1, &history);
150
151 exit(EXIT_SUCCESS);
152 }
153
154
155