1 /*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2021 Paul Ramsey <pramsey@cleverelephant.ca>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software Foundation,
20 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21 *
22 */
23
24 #include "../../postgis_config.h"
25 /* #define POSTGIS_DEBUG_LEVEL 4 */
26
27 #include "librtcore.h"
28 #include "librtcore_internal.h"
29
30 /******************************************************************************
31 * rt_raster_gdal_warp()
32 ******************************************************************************/
33
34 typedef struct
35 {
36 struct {
37 GDALDatasetH ds;
38 GDALDriverH drv;
39 int destroy_drv;
40 } src;
41
42 struct {
43 OGRSFDriverH drv;
44 OGRDataSourceH ds;
45 OGRLayerH lyr;
46 int srid;
47 OGRwkbGeometryType gtype;
48 } dst;
49
50 } _rti_contour_arg;
51
52
53 static void
_rti_contour_arg_init(_rti_contour_arg * arg)54 _rti_contour_arg_init(_rti_contour_arg* arg)
55 {
56 memset(arg, 0, sizeof(_rti_contour_arg));
57 return;
58 }
59
60
61 static int
_rti_contour_arg_destroy(_rti_contour_arg * arg)62 _rti_contour_arg_destroy(_rti_contour_arg* arg)
63 {
64 if(arg->src.ds != NULL)
65 GDALClose(arg->src.ds);
66
67 if (arg->src.drv != NULL && arg->src.destroy_drv) {
68 GDALDeregisterDriver(arg->src.drv);
69 GDALDestroyDriver(arg->src.drv);
70 }
71
72 if (arg->dst.ds != NULL)
73 OGR_DS_Destroy(arg->dst.ds);
74
75 return FALSE;
76 }
77
78
79
80 /**
81 * Return palloc'ed list of contours.
82 * @param src_raster : raster to generate contour from
83 * @param options : CSList of OPTION=VALUE strings for the
84 * contour routine, see https://gdal.org/api/gdal_alg.html?highlight=contour#_CPPv419GDALContourGenerate15GDALRasterBandHddiPdidPvii16GDALProgressFuncPv
85 * @param src_srs : Coordinate reference system string for raster
86 * @param ncontours : Output parameter for length of contour list
87 * @param contours : Output palloc'ed list of contours, caller to free
88 */
rt_raster_gdal_contour(rt_raster src_raster,int src_band,int src_srid,const char * src_srs,double contour_interval,double contour_base,int fixed_level_count,double * fixed_levels,int polygonize,size_t * ncontours,struct rt_contour_t ** contours)89 int rt_raster_gdal_contour(
90 /* input parameters */
91 rt_raster src_raster,
92 int src_band,
93 int src_srid,
94 const char* src_srs,
95 double contour_interval,
96 double contour_base,
97 int fixed_level_count,
98 double *fixed_levels,
99 int polygonize,
100 /* output parameters */
101 size_t *ncontours,
102 struct rt_contour_t **contours
103 )
104 {
105 CPLErr cplerr;
106 OGRErr ogrerr;
107 GDALRasterBandH hBand;
108 int nfeatures = 0, i = 0;
109 OGRFeatureH hFeat;
110
111 _rti_contour_arg arg;
112 _rti_contour_arg_init(&arg);
113
114 /* Load raster into GDAL memory */
115 arg.src.ds = rt_raster_to_gdal_mem(src_raster, src_srs, NULL, NULL, 0, &(arg.src.drv), &(arg.src.destroy_drv));
116 /* Extract the desired band */
117 hBand = GDALGetRasterBand(arg.src.ds, src_band);
118
119 /* Set up the OGR destination data store */
120 arg.dst.srid = src_srid;
121 arg.dst.drv = OGRGetDriverByName("Memory");
122 if (!arg.dst.drv)
123 return _rti_contour_arg_destroy(&arg);
124
125 arg.dst.ds = OGR_Dr_CreateDataSource(arg.dst.drv, "contour_ds", NULL);
126 if (!arg.dst.ds)
127 return _rti_contour_arg_destroy(&arg);
128
129 /* Polygonize is 2.4+ only */
130 #if POSTGIS_GDAL_VERSION >= 24
131 arg.dst.gtype = polygonize ? wkbPolygon : wkbLineString;
132 #else
133 arg.dst.gtype = wkbLineString;
134 #endif
135
136 /* Layer has geometry, elevation, id */
137 arg.dst.lyr = OGR_DS_CreateLayer(arg.dst.ds, "contours", NULL, arg.dst.gtype, NULL);
138 if (!arg.dst.lyr)
139 return _rti_contour_arg_destroy(&arg);
140
141 /* ID field */
142 OGRFieldDefnH hFldId = OGR_Fld_Create("id", OFTInteger);
143 ogrerr = OGR_L_CreateField(arg.dst.lyr, hFldId, TRUE);
144 if (ogrerr != OGRERR_NONE)
145 return _rti_contour_arg_destroy(&arg);
146
147 /* ELEVATION field */
148 OGRFieldDefnH hFldElevation = OGR_Fld_Create("elevation", OFTReal);
149 ogrerr = OGR_L_CreateField(arg.dst.lyr, hFldElevation, TRUE);
150 if (ogrerr != OGRERR_NONE)
151 return _rti_contour_arg_destroy(&arg);
152
153 // GDALContourGenerate( GDALRasterBandH hBand,
154 // double dfContourInterval, double dfContourBase,
155 // int nFixedLevelCount, double *padfFixedLevels,
156 // int bUseNoData, double dfNoDataValue,
157 // void *hLayer, int iIDField, int iElevField,
158 // GDALProgressFunc pfnProgress, void *pProgressArg );
159
160 int use_no_data = 0;
161 double no_data_value = GDALGetRasterNoDataValue(hBand, &use_no_data);;
162
163 /* Run the contouring routine, filling up the OGR layer */
164 cplerr = GDALContourGenerate(
165 hBand,
166 contour_interval, contour_base,
167 fixed_level_count, fixed_levels,
168 use_no_data, no_data_value,
169 arg.dst.lyr, 0, 1, // OGRLayer, idFieldNum, elevFieldNum
170 NULL, // GDALProgressFunc pfnProgress
171 NULL // void *pProgressArg
172 );
173
174 if (cplerr != CE_None)
175 return _rti_contour_arg_destroy(&arg);
176
177 /* Convert the OGR layer into PostGIS geometries */
178 nfeatures = OGR_L_GetFeatureCount(arg.dst.lyr, TRUE);
179 if (nfeatures < 0)
180 return _rti_contour_arg_destroy(&arg);
181
182 *contours = rtalloc(sizeof(struct rt_contour_t) * nfeatures);
183 OGR_L_ResetReading(arg.dst.lyr);
184 while ((hFeat = OGR_L_GetNextFeature(arg.dst.lyr))) {
185 size_t szWkb;
186 unsigned char *bufWkb;
187 struct rt_contour_t contour;
188 OGRGeometryH hGeom;
189 LWGEOM *geom;
190
191 /* Somehow we're still iterating, should not happen. */
192 if(i >= nfeatures) break;
193
194 /* Store elevation/id */
195 contour.id = OGR_F_GetFieldAsInteger(hFeat, 0);
196 contour.elevation = OGR_F_GetFieldAsDouble(hFeat, 1);
197 /* Convert OGR geometry to LWGEOM via WKB */
198 if (!(hGeom = OGR_F_GetGeometryRef(hFeat))) continue;
199 szWkb = OGR_G_WkbSize(hGeom);
200 bufWkb = rtalloc(szWkb);
201 if (OGR_G_ExportToWkb(hGeom, wkbNDR, bufWkb) != OGRERR_NONE) continue;
202 /* Reclaim feature and associated geometry memory */
203 OGR_F_Destroy(hFeat);
204 geom = lwgeom_from_wkb(bufWkb, szWkb, LW_PARSER_CHECK_NONE);
205 lwgeom_set_srid(geom, arg.dst.srid);
206 contour.geom = gserialized_from_lwgeom(geom, NULL);
207 lwgeom_free(geom);
208 rtdealloc(bufWkb);
209 (*contours)[i++] = contour;
210 }
211
212 /* Return total number of contours saved */
213 *ncontours = i;
214
215 /* Free all the non-database allocated structures */
216 _rti_contour_arg_destroy(&arg);
217
218 /* Done */
219 return TRUE;
220 }
221
222