window_map.c

Go to the documentation of this file.

#include <stdlib.h>
#include <grass/gis.h>
#include "G.h"


#define alloc_index(n) (COLUMN_MAPPING *) G_malloc((n)*sizeof(COLUMN_MAPPING))


int G__create_window_mapping(int fd)
{
    struct fileinfo *fcb = &G__.fileinfo[fd];
    COLUMN_MAPPING *col;
    int i;
    int x;
    double C1, C2;
    double west;

    G__init_window();

    if (fcb->open_mode >= 0 && fcb->open_mode != OPEN_OLD)      /* open for write? */
        return 0;
    if (fcb->open_mode == OPEN_OLD)     /* already open ? */
        G_free(fcb->col_map);

    col = fcb->col_map = alloc_index(G__.window.cols);

    /*
     * for each column in the window, go to center of the cell,
     * compute nearest column in the data file
     * if column is not in data file, set column to 0
     *
     * for lat/lon move window so that west is bigger than
     * cellhd west.
     */
    west = G__.window.west;
    if (G__.window.proj == PROJECTION_LL) {
        while (west > fcb->cellhd.west + 360.0)
            west -= 360.0;
        while (west < fcb->cellhd.west)
            west += 360.0;
    }

    C1 = G__.window.ew_res / fcb->cellhd.ew_res;
    C2 = (west - fcb->cellhd.west +
          G__.window.ew_res / 2.0) / fcb->cellhd.ew_res;
    for (i = 0; i < G__.window.cols; i++) {
        x = C2;
        if (C2 < x)             /* adjust for rounding of negatives */
            x--;
        if (x < 0 || x >= fcb->cellhd.cols)     /* not in data file */
            x = -1;
        *col++ = x + 1;
        C2 += C1;
    }

    /* do wrap around for lat/lon */
    if (G__.window.proj == PROJECTION_LL) {
        col = fcb->col_map;
        C2 = (west - 360.0 - fcb->cellhd.west +
              G__.window.ew_res / 2.0) / fcb->cellhd.ew_res;
        for (i = 0; i < G__.window.cols; i++) {
            x = C2;
            if (C2 < x)         /* adjust for rounding of negatives */
                x--;
            if (x < 0 || x >= fcb->cellhd.cols) /* not in data file */
                x = -1;
            if (*col == 0)      /* only change those not already set */
                *col = x + 1;
            col++;
            C2 += C1;
        }
    }

    G_debug(3, "create window mapping (%d columns)", G__.window.cols);
/*  for (i = 0; i < G__.window.cols; i++)
        fprintf(stderr, "%s%ld", i % 15 ? " " : "\n", (long)fcb->col_map[i]);
    fprintf(stderr, "\n");
*/

    /* compute C1,C2 for row window mapping */
    fcb->C1 = G__.window.ns_res / fcb->cellhd.ns_res;
    fcb->C2 =
        (fcb->cellhd.north - G__.window.north +
         G__.window.ns_res / 2.0) / fcb->cellhd.ns_res;

    return 0;
}


double G_northing_to_row(double north, const struct Cell_head *window)
{
    return (window->north - north) / window->ns_res;
}


double G_adjust_east_longitude(double east, double west)
{
    while (east > west + 360.0)
        east -= 360.0;
    while (east <= west)
        east += 360.0;

    return east;
}


double G_adjust_easting(double east, const struct Cell_head *window)
{
    if (window->proj == PROJECTION_LL) {
        east = G_adjust_east_longitude(east, window->west);
        if (east > window->east && east == window->west + 360)
            east = window->west;
    }

    return east;
}


double G_easting_to_col(double east, const struct Cell_head *window)
{
    east = G_adjust_easting(east, window);

    return (east - window->west) / window->ew_res;
}


double G_row_to_northing(double row, const struct Cell_head *window)
{
    return window->north - row * window->ns_res;
}


double G_col_to_easting(double col, const struct Cell_head *window)
{
    return window->west + col * window->ew_res;
}


int G_window_rows(void)
{
    G__init_window();

    return G__.window.rows;
}


int G_window_cols(void)
{
    G__init_window();

    return G__.window.cols;
}


int G__init_window(void)
{
    if (!G__.window_set) {
        G__.window_set = 1;
        G_get_window(&G__.window);
    }

    return 0;
}


int G_row_repeat_nomask(int fd, int row)
{
    struct fileinfo *fcb = &G__.fileinfo[fd];
    double f;
    int r1, r2;
    int count;

    count = 1;

    /* r1 is the row in the raster map itself.
     * r2 is the next row(s) in the raster map
     * see get_row.c for details on this calculation
     */
    f = row * fcb->C1 + fcb->C2;
    r1 = f;
    if (f < r1)
        r1--;

    while (++row < G__.window.rows) {
        f = row * fcb->C1 + fcb->C2;
        r2 = f;
        if (f < r2)
            r2--;
        if (r1 != r2)
            break;

        count++;
    }

    return count;
}