← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.proj.seg/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.proj.seg/main.c
1
 
 
2
 
/***************************************************************************
3
 
*
4
 
* MODULE:       r.proj
5
 
*
6
 
* AUTHOR(S):    Martin Schroeder
7
 
*                University of Heidelberg
8
 
*                Dept. of Geography
9
 
*                emes@geo0.geog.uni-heidelberg.de
10
 
*
11
 
*                (With the help of a lot of existing GRASS sources, in 
12
 
*                 particular v.proj) 
13
 
*
14
 
* PURPOSE:      r.proj converts a map to a new geographic projection. It reads a
15
 
*               map from a different location, projects it and write it out
16
 
*               to the current location. The projected data is resampled with
17
 
*               one of three different methods: nearest neighbor, bilinear and
18
 
*               cubic convolution.
19
 
*
20
 
* COPYRIGHT:    (C) 2001 by the GRASS Development Team
21
 
*
22
 
*               This program is free software under the GNU General Public
23
 
*               License (>=v2). Read the file COPYING that comes with GRASS
24
 
*               for details.
25
 
*
26
 
* Changes
27
 
*                Morten Hulden <morten@untamo.net>, Aug 2000:
28
 
*                - aborts if input map is outside current location.
29
 
*                - can handle projections (conic, azimuthal etc) where 
30
 
*                part of the map may fall into areas where south is 
31
 
*                upward and east is leftward.
32
 
*                - avoids passing location edge coordinates to PROJ
33
 
*                (they may be invalid in some projections).
34
 
*                - output map will be clipped to borders of the current region.
35
 
*                - output map cell edges and centers will coinside with those 
36
 
*                of the current region.
37
 
*                - output map resolution (unless changed explicitly) will
38
 
*                match (exactly) the resolution of the current region.
39
 
*                - if the input map is smaller than the current region, the 
40
 
*                output map will only cover the overlapping area.
41
 
*                - if the input map is larger than the current region, only the
42
 
*                needed amount of memory will be allocated for the projection
43
 
*       
44
 
*                Bugfixes 20050328: added floor() before (int) typecasts to in avoid
45
 
*                assymetrical rounding errors. Added missing offset outcellhd.ew_res/2 
46
 
*                to initial xcoord for each row in main projection loop (we want to  project 
47
 
*                center of cell, not border).
48
 
*
49
 
*                Glynn Clements 2006: Use G_interp_* functions, modified      
50
 
*                  version of r.proj which uses a tile cache instead  of loading the
51
 
*                  entire map into memory.
52
 
 
53
 
*****************************************************************************/
54
 
 
55
 
 
56
 
#include <stdio.h>
57
 
#include <stdlib.h>
58
 
#include <string.h>
59
 
#include <unistd.h>
60
 
#include <grass/gis.h>
61
 
#include <grass/gprojects.h>
62
 
#include <grass/glocale.h>
63
 
#include "r.proj.h"
64
 
 
65
 
/* modify this table to add new methods */
66
 
struct menu menu[] = {
67
 
    {p_nearest, "nearest", "nearest neighbor"},
68
 
    {p_bilinear, "bilinear", "bilinear"},
69
 
    {p_cubic, "cubic", "cubic convolution"},
70
 
    {NULL, NULL, NULL}
71
 
};
72
 
 
73
 
static char *make_ipol_list(void);
74
 
 
75
 
int main(int argc, char **argv)
76
 
{
77
 
    char *mapname,              /* ptr to name of output layer  */
78
 
     *setname,                  /* ptr to name of input mapset  */
79
 
     *ipolname;                 /* name of interpolation method */
80
 
 
81
 
    int fdi,                    /* input map file descriptor    */
82
 
      fdo,                      /* output map file descriptor   */
83
 
      method,                   /* position of method in table  */
84
 
      permissions,              /* mapset permissions           */
85
 
      cell_type,                /* output celltype              */
86
 
      cell_size,                /* size of a cell in bytes      */
87
 
      row, col,                 /* counters                     */
88
 
      irows, icols,             /* original rows, cols          */
89
 
      orows, ocols, have_colors,        /* Input map has a colour table */
90
 
      overwrite,                /* Overwrite                    */
91
 
      curr_proj;                /* output projection (see gis.h) */
92
 
 
93
 
    void *obuffer,              /* buffer that holds one output row     */
94
 
     *obufptr;                  /* column ptr in output buffer  */
95
 
    struct cache *ibuffer;      /* buffer that holds the input map      */
96
 
    func interpolate;           /* interpolation routine        */
97
 
 
98
 
    double xcoord1, xcoord2,    /* temporary x coordinates      */
99
 
      ycoord1, ycoord2,         /* temporary y coordinates      */
100
 
      col_idx,                  /* column index in input matrix */
101
 
      row_idx,                  /* row index in input matrix    */
102
 
      onorth, osouth,           /* save original border coords  */
103
 
      oeast, owest, inorth, isouth, ieast, iwest;
104
 
    char north_str[30], south_str[30], east_str[30], west_str[30];
105
 
 
106
 
    struct Colors colr;         /* Input map colour table       */
107
 
    struct History history;
108
 
 
109
 
    struct pj_info iproj,       /* input map proj parameters    */
110
 
      oproj;                    /* output map proj parameters   */
111
 
 
112
 
    struct Key_Value *in_proj_info,     /* projection information of    */
113
 
     *in_unit_info,             /* input and output mapsets     */
114
 
     *out_proj_info, *out_unit_info;
115
 
 
116
 
    struct GModule *module;
117
 
 
118
 
    struct Flag *list,          /* list files in source location */
119
 
     *nocrop,                   /* don't crop output map        */
120
 
     *print_bounds,             /* print output bounds and exit */
121
 
     *gprint_bounds;            /* same but print shell style   */
122
 
 
123
 
    struct Option *imapset,     /* name of input mapset         */
124
 
     *inmap,                    /* name of input layer          */
125
 
     *inlocation,               /* name of input location       */
126
 
     *outmap,                   /* name of output layer         */
127
 
     *indbase,                  /* name of input database       */
128
 
     *interpol,                 /* interpolation method:
129
 
                                   nearest neighbor, bilinear, cubic */
130
 
     *memory,                   /* amount of memory for cache   */
131
 
     *res;                      /* resolution of target map     */
132
 
    struct Cell_head incellhd,  /* cell header of input map     */
133
 
      outcellhd;                /* and output map               */
134
 
 
135
 
 
136
 
    G_gisinit(argv[0]);
137
 
 
138
 
    module = G_define_module();
139
 
    module->keywords = _("raster, projection, transformation");
140
 
     module->description =
141
 
        _("Re-projects a raster map from one location to the current location.");
142
 
 
143
 
    inmap = G_define_standard_option(G_OPT_R_INPUT);
144
 
    inmap->description = _("Name of input raster map to re-project");
145
 
    inmap->required = NO;
146
 
    inmap->guisection = _("Source");
147
 
    
148
 
    inlocation = G_define_option();
149
 
    inlocation->key = "location";
150
 
    inlocation->type = TYPE_STRING;
151
 
    inlocation->required = YES;
152
 
    inlocation->description = _("Location containing input raster map");
153
 
    inlocation->gisprompt = "old,location,location";
154
 
    inlocation->key_desc = "name";
155
 
 
156
 
    imapset = G_define_option();
157
 
    imapset->key = "mapset";
158
 
    imapset->type = TYPE_STRING;
159
 
    imapset->required = NO;
160
 
    imapset->description = _("Mapset containing input raster map");
161
 
    imapset->gisprompt = "old,mapset,mapset";
162
 
    imapset->key_desc = "name";
163
 
    imapset->guisection = _("Source");
164
 
 
165
 
    indbase = G_define_option();
166
 
    indbase->key = "dbase";
167
 
    indbase->type = TYPE_STRING;
168
 
    indbase->required = NO;
169
 
    indbase->description = _("Path to GRASS database of input location");
170
 
    indbase->gisprompt = "old,dbase,dbase";
171
 
    indbase->key_desc = "path";
172
 
    indbase->guisection = _("Source");
173
 
 
174
 
    outmap = G_define_standard_option(G_OPT_R_OUTPUT);
175
 
    outmap->required = NO;
176
 
    outmap->description = _("Name for output raster map (default: input)");
177
 
    outmap->guisection = _("Target");
178
 
    
179
 
    ipolname = make_ipol_list();
180
 
 
181
 
    interpol = G_define_option();
182
 
    interpol->key = "method";
183
 
    interpol->type = TYPE_STRING;
184
 
    interpol->required = NO;
185
 
    interpol->answer = "nearest";
186
 
    interpol->options = ipolname;
187
 
    interpol->description = _("Interpolation method to use");
188
 
    interpol->guisection = _("Target");
189
 
 
190
 
    memory = G_define_option();
191
 
    memory->key = "memory";
192
 
    memory->type = TYPE_INTEGER;
193
 
    memory->required = NO;
194
 
    memory->description = _("Cache size (MiB)");
195
 
 
196
 
    res = G_define_option();
197
 
    res->key = "resolution";
198
 
    res->type = TYPE_DOUBLE;
199
 
    res->required = NO;
200
 
    res->description = _("Resolution of output map");
201
 
    res->guisection = _("Target");
202
 
 
203
 
    list = G_define_flag();
204
 
    list->key = 'l';
205
 
    list->description = _("List raster maps in input location and exit");
206
 
 
207
 
    nocrop = G_define_flag();
208
 
    nocrop->key = 'n';
209
 
    nocrop->description = _("Do not perform region cropping optimization");
210
 
 
211
 
    print_bounds = G_define_flag();
212
 
    print_bounds->key = 'p';
213
 
    print_bounds->description =
214
 
        _("Print input map's bounds in the current projection and exit");
215
 
    print_bounds->guisection = _("Target");
216
 
    
217
 
    gprint_bounds = G_define_flag();
218
 
    gprint_bounds->key = 'g';
219
 
    gprint_bounds->description =
220
 
        _("Print input map's bounds in the current projection and exit (shell style)");
221
 
    gprint_bounds->guisection = _("Target");
222
 
 
223
 
    /* The parser checks if the map already exists in current mapset,
224
 
       we switch out the check and do it
225
 
       in the module after the parser */
226
 
    overwrite = G_check_overwrite(argc, argv);
227
 
 
228
 
    if (G_parser(argc, argv))
229
 
        exit(EXIT_FAILURE);
230
 
 
231
 
 
232
 
    /* get the method */
233
 
    for (method = 0; (ipolname = menu[method].name); method++)
234
 
        if (strcmp(ipolname, interpol->answer) == 0)
235
 
            break;
236
 
 
237
 
    if (!ipolname)
238
 
        G_fatal_error(_("<%s=%s> unknown %s"),
239
 
                      interpol->key, interpol->answer, interpol->key);
240
 
    interpolate = menu[method].method;
241
 
 
242
 
    mapname = outmap->answer ? outmap->answer : inmap->answer;
243
 
    if (mapname && !list->answer && !overwrite &&
244
 
        G_find_cell(mapname, G_mapset()))
245
 
        G_fatal_error(_("option <%s>: <%s> exists."), "output", mapname);
246
 
 
247
 
    setname = imapset->answer ? imapset->answer : G_store(G_mapset());
248
 
    if (strcmp(inlocation->answer, G_location()) == 0 &&
249
 
        (!indbase->answer || strcmp(indbase->answer, G_gisdbase()) == 0))
250
 
#if 0
251
 
        G_fatal_error(_("Input and output locations can not be the same"));
252
 
#else
253
 
        G_warning(_("Input and output locations are the same"));
254
 
#endif
255
 
    G_get_window(&outcellhd);
256
 
 
257
 
    if(gprint_bounds->answer && !print_bounds->answer)
258
 
        print_bounds->answer = gprint_bounds->answer;
259
 
    curr_proj = G_projection();
260
 
 
261
 
    /* Get projection info for output mapset */
262
 
    if ((out_proj_info = G_get_projinfo()) == NULL)
263
 
        G_fatal_error(_("Unable to get projection info of output raster map"));
264
 
 
265
 
    if ((out_unit_info = G_get_projunits()) == NULL)
266
 
        G_fatal_error(_("Unable to get projection units of output raster map"));
267
 
 
268
 
    if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
269
 
        G_fatal_error(_("Unable to get projection key values of output raster map"));
270
 
 
271
 
    /* Change the location           */
272
 
    G__create_alt_env();
273
 
    G__setenv("GISDBASE", indbase->answer ? indbase->answer : G_gisdbase());
274
 
    G__setenv("LOCATION_NAME", inlocation->answer);
275
 
 
276
 
    permissions = G__mapset_permissions(setname);
277
 
    if (permissions < 0)        /* can't access mapset       */
278
 
        G_fatal_error(_("Mapset <%s> in input location <%s> - %s"),
279
 
                      setname, inlocation->answer,
280
 
                      permissions == 0 ? _("permission denied")
281
 
                      : _("not found"));
282
 
 
283
 
    /* if requested, list the raster maps in source location - MN 5/2001 */
284
 
    if (list->answer) {
285
 
        int i;
286
 
        char **list;
287
 
        G_verbose_message(_("Checking location <%s> mapset <%s>"),
288
 
                          inlocation->answer, setname);
289
 
        list = G_list(G_ELEMENT_RASTER, G__getenv("GISDBASE"),
290
 
                      G__getenv("LOCATION_NAME"), setname);
291
 
        for (i = 0; list[i]; i++) {
292
 
            fprintf(stdout, "%s\n", list[i]);
293
 
        }
294
 
        fflush(stdout);
295
 
        exit(EXIT_SUCCESS);     /* leave v.proj after listing */
296
 
    }
297
 
 
298
 
    if (!inmap->answer)
299
 
        G_fatal_error(_("Required parameter <%s> not set"), inmap->key);
300
 
 
301
 
    if (!G_find_cell(inmap->answer, setname))
302
 
        G_fatal_error(_("Raster map <%s> in location <%s> in mapset <%s> not found"),
303
 
                      inmap->answer, inlocation->answer, setname);
304
 
 
305
 
    /* Read input map colour table */
306
 
    have_colors = G_read_colors(inmap->answer, setname, &colr);
307
 
 
308
 
    /* Get projection info for input mapset */
309
 
    if ((in_proj_info = G_get_projinfo()) == NULL)
310
 
        G_fatal_error(_("Unable to get projection info of input map"));
311
 
 
312
 
    if ((in_unit_info = G_get_projunits()) == NULL)
313
 
        G_fatal_error(_("Unable to get projection units of input map"));
314
 
 
315
 
    if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
316
 
        G_fatal_error(_("Unable to get projection key values of input map"));
317
 
 
318
 
    G_free_key_value(in_proj_info);
319
 
    G_free_key_value(in_unit_info);
320
 
    G_free_key_value(out_proj_info);
321
 
    G_free_key_value(out_unit_info);
322
 
    pj_print_proj_params(&iproj, &oproj);
323
 
 
324
 
    /* this call causes r.proj to read the entire map into memeory */
325
 
    G_get_cellhd(inmap->answer, setname, &incellhd);
326
 
 
327
 
    G_set_window(&incellhd);
328
 
 
329
 
    if (G_projection() == PROJECTION_XY)
330
 
        G_fatal_error(_("Unable to work with unprojected data (xy location)"));
331
 
 
332
 
    /* Save default borders so we can show them later */
333
 
    inorth = incellhd.north;
334
 
    isouth = incellhd.south;
335
 
    ieast = incellhd.east;
336
 
    iwest = incellhd.west;
337
 
    irows = incellhd.rows;
338
 
    icols = incellhd.cols;
339
 
 
340
 
    onorth = outcellhd.north;
341
 
    osouth = outcellhd.south;
342
 
    oeast = outcellhd.east;
343
 
    owest = outcellhd.west;
344
 
    orows = outcellhd.rows;
345
 
    ocols = outcellhd.cols;
346
 
 
347
 
 
348
 
    if (print_bounds->answer) {
349
 
        G_message(_("Input map <%s@%s> in location <%s>:"),
350
 
            inmap->answer, setname, inlocation->answer);
351
 
 
352
 
        if (pj_do_proj(&iwest, &isouth, &iproj, &oproj) < 0)
353
 
            G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
354
 
        if (pj_do_proj(&ieast, &inorth, &iproj, &oproj) < 0)
355
 
            G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
356
 
 
357
 
        G_format_northing(inorth, north_str, curr_proj);
358
 
        G_format_northing(isouth, south_str, curr_proj);
359
 
        G_format_easting(ieast, east_str, curr_proj);
360
 
        G_format_easting(iwest, west_str, curr_proj);
361
 
 
362
 
        if(gprint_bounds->answer) {
363
 
            fprintf(stdout, "n=%s s=%s w=%s e=%s rows=%d cols=%d\n",
364
 
                north_str, south_str, west_str, east_str, irows, icols);
365
 
        }
366
 
        else {
367
 
            fprintf(stdout, "Source cols: %d\n", icols);
368
 
            fprintf(stdout, "Source rows: %d\n", irows);
369
 
            fprintf(stdout, "Local north: %s\n",  north_str);
370
 
            fprintf(stdout, "Local south: %s\n", south_str);
371
 
            fprintf(stdout, "Local west: %s\n", west_str);
372
 
            fprintf(stdout, "Local east: %s\n", east_str);
373
 
        }
374
 
 
375
 
        /* somehow approximate local ewres, nsres ?? (use 'g.region -m' on lat/lon side) */
376
 
 
377
 
        exit(EXIT_SUCCESS);
378
 
    }
379
 
 
380
 
 
381
 
    /* Cut non-overlapping parts of input map */
382
 
    if (!nocrop->answer)
383
 
        bordwalk(&outcellhd, &incellhd, &oproj, &iproj);
384
 
 
385
 
    /* Add 2 cells on each side for bilinear/cubic & future interpolation methods */
386
 
    /* (should probably be a factor based on input and output resolution) */
387
 
    incellhd.north += 2 * incellhd.ns_res;
388
 
    incellhd.east += 2 * incellhd.ew_res;
389
 
    incellhd.south -= 2 * incellhd.ns_res;
390
 
    incellhd.west -= 2 * incellhd.ew_res;
391
 
    if (incellhd.north > inorth)
392
 
        incellhd.north = inorth;
393
 
    if (incellhd.east > ieast)
394
 
        incellhd.east = ieast;
395
 
    if (incellhd.south < isouth)
396
 
        incellhd.south = isouth;
397
 
    if (incellhd.west < iwest)
398
 
        incellhd.west = iwest;
399
 
 
400
 
    G_set_window(&incellhd);
401
 
 
402
 
    /* And switch back to original location */
403
 
 
404
 
    G__switch_env();
405
 
 
406
 
    /* Adjust borders of output map */
407
 
 
408
 
    if (!nocrop->answer)
409
 
        bordwalk(&incellhd, &outcellhd, &iproj, &oproj);
410
 
 
411
 
#if 0
412
 
    outcellhd.west = outcellhd.south = HUGE_VAL;
413
 
    outcellhd.east = outcellhd.north = -HUGE_VAL;
414
 
    for (row = 0; row < incellhd.rows; row++) {
415
 
        ycoord1 = G_row_to_northing((double)(row + 0.5), &incellhd);
416
 
        for (col = 0; col < incellhd.cols; col++) {
417
 
            xcoord1 = G_col_to_easting((double)(col + 0.5), &incellhd);
418
 
            pj_do_proj(&xcoord1, &ycoord1, &iproj, &oproj);
419
 
            if (xcoord1 > outcellhd.east)
420
 
                outcellhd.east = xcoord1;
421
 
            if (ycoord1 > outcellhd.north)
422
 
                outcellhd.north = ycoord1;
423
 
            if (xcoord1 < outcellhd.west)
424
 
                outcellhd.west = xcoord1;
425
 
            if (ycoord1 < outcellhd.south)
426
 
                outcellhd.south = ycoord1;
427
 
        }
428
 
    }
429
 
#endif
430
 
 
431
 
    if (res->answer != NULL)    /* set user defined resolution */
432
 
        outcellhd.ns_res = outcellhd.ew_res = atof(res->answer);
433
 
 
434
 
    G_adjust_Cell_head(&outcellhd, 0, 0);
435
 
    G_set_window(&outcellhd);
436
 
 
437
 
    G_message(" ");
438
 
    G_message(_("Input:"));
439
 
    G_message(_("Cols: %d (%d)"), incellhd.cols, icols);
440
 
    G_message(_("Rows: %d (%d)"), incellhd.rows, irows);
441
 
    G_message(_("North: %f (%f)"), incellhd.north, inorth);
442
 
    G_message(_("South: %f (%f)"), incellhd.south, isouth);
443
 
    G_message(_("West: %f (%f)"), incellhd.west, iwest);
444
 
    G_message(_("East: %f (%f)"), incellhd.east, ieast);
445
 
    G_message(_("EW-res: %f"), incellhd.ew_res);
446
 
    G_message(_("NS-res: %f"), incellhd.ns_res);
447
 
 
448
 
    G_message(" ");
449
 
    G_message(_("Output:"));
450
 
    G_message(_("Cols: %d (%d)"), outcellhd.cols, ocols);
451
 
    G_message(_("Rows: %d (%d)"), outcellhd.rows, orows);
452
 
    G_message(_("North: %f (%f)"), outcellhd.north, onorth);
453
 
    G_message(_("South: %f (%f)"), outcellhd.south, osouth);
454
 
    G_message(_("West: %f (%f)"), outcellhd.west, owest);
455
 
    G_message(_("East: %f (%f)"), outcellhd.east, oeast);
456
 
    G_message(_("EW-res: %f"), outcellhd.ew_res);
457
 
    G_message(_("NS-res: %f"), outcellhd.ns_res);
458
 
    G_message(" ");
459
 
 
460
 
    /* open and read the relevant parts of the input map and close it */
461
 
    G__switch_env();
462
 
    G_set_window(&incellhd);
463
 
    fdi = G_open_cell_old(inmap->answer, setname);
464
 
    cell_type = G_get_raster_map_type(fdi);
465
 
    ibuffer = readcell(fdi, memory->answer);
466
 
    G_close_cell(fdi);
467
 
 
468
 
    G__switch_env();
469
 
    G_set_window(&outcellhd);
470
 
 
471
 
    if (strcmp(interpol->answer, "nearest") == 0) {
472
 
        fdo = G_open_raster_new(mapname, cell_type);
473
 
        obuffer = (CELL *) G_allocate_raster_buf(cell_type);
474
 
    }
475
 
    else {
476
 
        fdo = G_open_fp_cell_new(mapname);
477
 
        cell_type = FCELL_TYPE;
478
 
        obuffer = (FCELL *) G_allocate_raster_buf(cell_type);
479
 
    }
480
 
 
481
 
    cell_size = G_raster_size(cell_type);
482
 
 
483
 
    xcoord1 = xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
484
 
     /**/ ycoord1 = ycoord2 = outcellhd.north - (outcellhd.ns_res / 2);
485
 
     /**/ G_important_message(_("Projecting..."));
486
 
    G_percent(0, outcellhd.rows, 2);
487
 
 
488
 
    for (row = 0; row < outcellhd.rows; row++) {
489
 
        obufptr = obuffer;
490
 
 
491
 
        for (col = 0; col < outcellhd.cols; col++) {
492
 
            /* project coordinates in output matrix to       */
493
 
            /* coordinates in input matrix                   */
494
 
            if (pj_do_proj(&xcoord1, &ycoord1, &oproj, &iproj) < 0)
495
 
                G_set_null_value(obufptr, 1, cell_type);
496
 
            else {
497
 
                /* convert to row/column indices of input matrix */
498
 
                col_idx = (xcoord1 - incellhd.west) / incellhd.ew_res;
499
 
                row_idx = (incellhd.north - ycoord1) / incellhd.ns_res;
500
 
 
501
 
                /* and resample data point               */
502
 
                interpolate(ibuffer, obufptr, cell_type,
503
 
                            &col_idx, &row_idx, &incellhd);
504
 
            }
505
 
 
506
 
            obufptr = G_incr_void_ptr(obufptr, cell_size);
507
 
            xcoord2 += outcellhd.ew_res;
508
 
            xcoord1 = xcoord2;
509
 
            ycoord1 = ycoord2;
510
 
        }
511
 
 
512
 
        if (G_put_raster_row(fdo, obuffer, cell_type) < 0)
513
 
            G_fatal_error(_("Failed writing raster map <%s> row %d"), mapname,
514
 
                          row);
515
 
 
516
 
        xcoord1 = xcoord2 = outcellhd.west + (outcellhd.ew_res / 2);
517
 
        ycoord2 -= outcellhd.ns_res;
518
 
        ycoord1 = ycoord2;
519
 
        G_percent(row, outcellhd.rows - 1, 2);
520
 
    }
521
 
 
522
 
    G_close_cell(fdo);
523
 
 
524
 
    if (have_colors > 0) {
525
 
        G_write_colors(mapname, G_mapset(), &colr);
526
 
        G_free_colors(&colr);
527
 
    }
528
 
 
529
 
    G_short_history(mapname, "raster", &history);
530
 
    G_command_history(&history);
531
 
    G_write_history(mapname, &history);
532
 
 
533
 
    G_done_msg("");
534
 
    exit(EXIT_SUCCESS);
535
 
}
536
 
 
537
 
static char *make_ipol_list(void)
538
 
{
539
 
    int size = 0;
540
 
    int i;
541
 
    char *buf;
542
 
 
543
 
    for (i = 0; menu[i].name; i++)
544
 
        size += strlen(menu[i].name) + 1;
545
 
 
546
 
    buf = G_malloc(size);
547
 
    *buf = '\0';
548
 
 
549
 
    for (i = 0; menu[i].name; i++) {
550
 
        if (i)
551
 
            strcat(buf, ",");
552
 
        strcat(buf, menu[i].name);
553
 
    }
554
 
 
555
 
    return buf;
556
 
}