melfp/outlets.c at main · HydroCSLab/melfp · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include "global.h"

struct outlet_list *read_outlets(const char *outlets_path,
                                 const char *layer_name, const char *opts,
                                 const char *id_col,
                                 struct raster_map *dir_map, int find_full)
{
    struct outlet_list *outlet_l = malloc(sizeof *outlet_l);
    const char **ds_opts = NULL;
    GDALDatasetH dataset;
    OGRLayerH layer;
    OGRFeatureDefnH feature_def;
    int id_field_idx;
    OGRFeatureH feature;

    if (opts) {
        char *o, *p;
        int n;

        strcpy((o = malloc(strlen(opts) + 1)), opts);
        for (p = o, n = 2; *p; p++)
            if (*p == ',')
                n++;
        ds_opts = malloc(sizeof *ds_opts * n);

        for (p = o, n = 0; *p; p++) {
            if (p == o || *(p - 1) == ',') {
                ds_opts[n] = p;
                if (p > o)
                    *(p - 1) = 0;
            }
        }
    }

    if (!(dataset =
          GDALOpenEx(outlets_path, GDAL_OF_VECTOR, NULL, ds_opts, NULL)))
        return NULL;

    if (!GDALDatasetGetLayerCount(dataset))
        return NULL;

    if (!(layer =
          layer_name ? GDALDatasetGetLayerByName(dataset, layer_name) :
          GDALDatasetGetLayer(dataset, 0)))
        return NULL;

    init_outlet_list(outlet_l);

    feature_def = OGR_L_GetLayerDefn(layer);
    if (strcmp(id_col, OGR_L_GetFIDColumn(layer)) == 0)
        id_field_idx = -1;
    else {
        OGRFieldDefnH id_field_def;

        if ((id_field_idx = OGR_FD_GetFieldIndex(feature_def, id_col)) < 0)
            return NULL;

        id_field_def = OGR_FD_GetFieldDefn(feature_def, id_field_idx);
        if (OGR_Fld_GetType(id_field_def) != OFTInteger)
            return NULL;
    }

    while ((feature = OGR_L_GetNextFeature(layer))) {
        OGRGeometryH geometry;
        int id;
        double x, y;

        geometry = OGR_F_GetGeometryRef(feature);
        if (geometry &&
            wkbFlatten(OGR_G_GetGeometryType(geometry)) == wkbPoint) {
            int row, col;

            x = OGR_G_GetX(geometry, 0);
            y = OGR_G_GetY(geometry, 0);
            id = id_field_idx >= 0 ?
                OGR_F_GetFieldAsInteger(feature, id_field_idx) :
                OGR_F_GetFID(feature);

            calc_row_col(dir_map, x, y, &row, &col);

            /* if the outlet is outside the computational region, skip */
            if (row >= 0 && row < dir_map->nrows && col >= 0 &&
                col < dir_map->ncols)
                add_outlet(outlet_l, row, col, id, find_full);
            else
                printf
                    ("Skip outlet %d at (%d, %d) outside the current region\n",
                     id, row, col);
        }

        OGR_F_Destroy(feature);
    }

    GDALClose(dataset);

    return outlet_l;
}

int write_outlets(const char *outlets_path, struct outlet_list *outlet_l)
{
    FILE *fp;
    int i;

    if (!(fp = fopen(outlets_path, "w")))
        return 1;

    for (i = 0; i < outlet_l->n; i++)
        fprintf(fp, "%d %d %d\n", outlet_l->row[i], outlet_l->col[i],
                outlet_l->id[i]);

    fclose(fp);

    return 0;
}