Skip to content

Commit 29e4ae8

Browse files
geographikageographika
authored
Merge pull request MapServer#7343 from geographika/raster-query-double
Use doubles for raster querying instead of float for Int32 support
1 parent b7fb8e4 commit 29e4ae8

File tree

6 files changed

+133
-19
lines changed

6 files changed

+133
-19
lines changed

‎msautotest/wxs/data/float64.tif‎

20.8 KB
Binary file not shown.

‎msautotest/wxs/data/int32.tif‎

10.6 KB
Binary file not shown.

‎msautotest/wxs/expected/wms_getfeatureinfo_float64_raster.json‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
{
2+
"type": "FeatureCollection",
3+
"name": "float64",
4+
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2157" } },
5+
"features": [
6+
{ "type": "Feature", "properties": { "x": "600182.31", "y": "628802.89", "value_0": 114140618.0, "value_list": "114140618", "red": "255", "green": "255", "blue": "255" }, "geometry": { "type": "Point", "coordinates": [ 600182.310000000055879, 628802.89000000001397 ] } }
7+
]
8+
}

‎msautotest/wxs/expected/wms_getfeatureinfo_int32_raster.json‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
{
2+
"type": "FeatureCollection",
3+
"name": "int32",
4+
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2157" } },
5+
"features": [
6+
{ "type": "Feature", "properties": { "x": "600182.31", "y": "628802.89", "value_0": 114140618, "value_list": "114140618", "red": "255", "green": "255", "blue": "255" }, "geometry": { "type": "Point", "coordinates": [ 600182.310000000055879, 628802.89000000001397 ] } }
7+
]
8+
}

‎msautotest/wxs/wms_getfeatureinfo_datatypes_raster.map‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#
2+
# Test WMS GetfeatureInfo with a raster Int32 layer
3+
#
4+
# REQUIRES: SUPPORTS=WMS
5+
#
6+
# RUN_PARMS: wms_getfeatureinfo_int32_raster.json [MAPSERV] "QUERY_STRING=map=[MAPFILE]&service=WMS&request=GetFeatureInfo&version=1.3.0&CRS=EPSG:2157&width=200&height=200&layers=int32&bbox=600181.89,628797.70,600191.27,628807.09&query_layers=int32&i=50&j=50&STYLES=&info_format=geojson" > [RESULT_DEMIME]
7+
# RUN_PARMS: wms_getfeatureinfo_float64_raster.json [MAPSERV] "QUERY_STRING=map=[MAPFILE]&service=WMS&request=GetFeatureInfo&version=1.3.0&CRS=EPSG:2157&width=200&height=200&layers=float64&bbox=600181.89,628797.70,600191.27,628807.09&query_layers=float64&i=50&j=50&STYLES=&info_format=geojson" > [RESULT_DEMIME]
8+
9+
MAP
10+
11+
NAME 'test'
12+
EXTENT 600029.81 628680.39 600319.81 628905.39
13+
SIZE 600 600
14+
IMAGETYPE PNG
15+
PROJECTION
16+
"init=epsg:2157"
17+
END
18+
19+
WEB
20+
METADATA
21+
"ows_enable_request" "*"
22+
"ows_srs" "EPSG:2157"
23+
"ows_title" "WMS Getfeatureinfo"
24+
"wms_getfeatureinfo_formatlist" "geojson"
25+
"wms_onlineresource" "http://localhost/cgi-bin/mapserv?map=mymap.map"
26+
END
27+
END
28+
29+
OUTPUTFORMAT
30+
NAME "geojson"
31+
DRIVER "OGR/GEOJSON"
32+
MIMETYPE "application/json; subtype=geojson; charset=utf-8"
33+
FORMATOPTION "FORM=SIMPLE"
34+
FORMATOPTION "STORAGE=memory"
35+
FORMATOPTION "LCO:NATIVE_MEDIA_TYPE=application/vnd.geo+json"
36+
FORMATOPTION "USE_FEATUREID=true"
37+
END
38+
39+
40+
LAYER
41+
NAME "int32"
42+
STATUS OFF
43+
TYPE RASTER
44+
TEMPLATE "blank"
45+
DATA "data/int32.tif"
46+
PROJECTION
47+
"EPSG:2157"
48+
END
49+
METADATA
50+
"gml_include_items" "all"
51+
"gml_value_0_type" "Long"
52+
END
53+
END
54+
55+
LAYER
56+
NAME "float64"
57+
STATUS OFF
58+
TYPE RASTER
59+
TEMPLATE "blank"
60+
DATA "data/float64.tif"
61+
PROJECTION
62+
"EPSG:2157"
63+
END
64+
METADATA
65+
"gml_include_items" "all"
66+
"gml_value_0_type" "Real"
67+
END
68+
END
69+
70+
END # of map file

‎src/maprasterquery.c‎

Lines changed: 47 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ typedef struct {
5959
double *qc_y;
6060
double *qc_x_reproj;
6161
double *qc_y_reproj;
62-
float *qc_values;
62+
double *qc_values;
6363
int *qc_class;
6464
int *qc_red;
6565
int *qc_green;
@@ -76,7 +76,7 @@ typedef struct {
7676
pointObj target_point;
7777

7878
GDALColorTableH hCT;
79-
79+
GDALDataType eDataType;
8080
double shape_tolerance;
8181

8282
} rasterLayerInfo;
@@ -211,7 +211,7 @@ static void msRasterLayerInfoInitialize(layerObj *layer)
211211
/************************************************************************/
212212

213213
static void msRasterQueryAddPixel(layerObj *layer, pointObj *location,
214-
pointObj *reprojectedLocation, float *values)
214+
pointObj *reprojectedLocation, double *values)
215215

216216
{
217217
rasterLayerInfo *rlinfo = (rasterLayerInfo *)layer->layerinfo;
@@ -238,8 +238,8 @@ static void msRasterQueryAddPixel(layerObj *layer, pointObj *location,
238238
(double *)msSmallCalloc(sizeof(double), rlinfo->query_alloc_max);
239239
rlinfo->qc_y_reproj =
240240
(double *)msSmallCalloc(sizeof(double), rlinfo->query_alloc_max);
241-
rlinfo->qc_values = (float *)msSmallCalloc(
242-
sizeof(float),
241+
rlinfo->qc_values = (double *)msSmallCalloc(
242+
sizeof(double),
243243
((size_t)rlinfo->query_alloc_max) * rlinfo->band_count);
244244
rlinfo->qc_red =
245245
(int *)msSmallCalloc(sizeof(int), rlinfo->query_alloc_max);
@@ -285,7 +285,7 @@ static void msRasterQueryAddPixel(layerObj *layer, pointObj *location,
285285
if (rlinfo->qc_values != NULL)
286286
rlinfo->qc_values = msSmallRealloc(
287287
rlinfo->qc_values,
288-
sizeof(float) * rlinfo->query_alloc_max * rlinfo->band_count);
288+
sizeof(double) * rlinfo->query_alloc_max * rlinfo->band_count);
289289
if (rlinfo->qc_class != NULL)
290290
rlinfo->qc_class = msSmallRealloc(rlinfo->qc_class,
291291
sizeof(int) * rlinfo->query_alloc_max);
@@ -385,7 +385,7 @@ static void msRasterQueryAddPixel(layerObj *layer, pointObj *location,
385385
/* -------------------------------------------------------------------- */
386386
if (rlinfo->qc_values != NULL)
387387
memcpy(rlinfo->qc_values + rlinfo->query_results * rlinfo->band_count,
388-
values, sizeof(float) * rlinfo->band_count);
388+
values, sizeof(double) * rlinfo->band_count);
389389

390390
/* -------------------------------------------------------------------- */
391391
/* Add to the results cache. */
@@ -408,7 +408,7 @@ static int msRasterQueryByRectLow(mapObj *map, layerObj *layer,
408408
double dfXMin, dfYMin, dfXMax, dfYMax, dfX, dfY, dfAdjustedRange;
409409
int nWinXOff, nWinYOff, nWinXSize, nWinYSize;
410410
int nRXSize, nRYSize;
411-
float *pafRaster;
411+
double *pafRaster;
412412
int nBandCount, *panBandMap, iPixel, iLine;
413413
CPLErr eErr;
414414
rasterLayerInfo *rlinfo;
@@ -504,15 +504,20 @@ static int msRasterQueryByRectLow(mapObj *map, layerObj *layer,
504504
/* band in the file. Later we will deal with the various band */
505505
/* selection criteria. */
506506
/* -------------------------------------------------------------------- */
507-
pafRaster = (float *)calloc(((size_t)nWinXSize) * nWinYSize * nBandCount,
508-
sizeof(float));
509-
MS_CHECK_ALLOC(pafRaster, sizeof(float) * nWinXSize * nWinYSize * nBandCount,
510-
-1);
511507

512-
eErr = GDALDatasetRasterIO(hDS, GF_Read, nWinXOff, nWinYOff, nWinXSize,
513-
nWinYSize, pafRaster, nWinXSize, nWinYSize,
514-
GDT_Float32, nBandCount, panBandMap,
515-
4 * nBandCount, 4 * nBandCount * nWinXSize, 4);
508+
size_t nPixels = (size_t)nWinXSize * nWinYSize * nBandCount;
509+
pafRaster = (double *)calloc(nPixels, sizeof(double));
510+
MS_CHECK_ALLOC(pafRaster, sizeof(double) * nPixels, -1);
511+
512+
// read raster block as GDT_Float64
513+
eErr = GDALDatasetRasterIO(
514+
hDS, GF_Read, nWinXOff, nWinYOff, nWinXSize, nWinYSize, pafRaster,
515+
nWinXSize, nWinYSize, GDT_Float64, nBandCount, panBandMap,
516+
sizeof(double) * nBandCount, sizeof(double) * nBandCount * nWinXSize,
517+
sizeof(double));
518+
519+
// store the datatype of the original raster to use for output formatting
520+
rlinfo->eDataType = GDALGetRasterDataType(GDALGetRasterBand(hDS, 1));
516521

517522
if (eErr != CE_None) {
518523
msSetError(MS_IOERR, "GDALDatasetRasterIO() failed: %s",
@@ -1179,13 +1184,36 @@ int msRASTERLayerGetShape(layerObj *layer, shapeObj *shape, resultObj *record) {
11791184
if (iValue != 0)
11801185
strlcat(szWork, ",", bufferSize);
11811186

1182-
snprintf(szWork + strlen(szWork), bufferSize - strlen(szWork), "%.8g",
1187+
snprintf(szWork + strlen(szWork), bufferSize - strlen(szWork),
1188+
"%.10g",
11831189
rlinfo->qc_values[shapeindex * rlinfo->band_count + iValue]);
11841190
}
11851191
} else if (EQUALN(layer->items[i], "value_", 6) && rlinfo->qc_values) {
11861192
int iValue = atoi(layer->items[i] + 6);
1187-
snprintf(szWork, bufferSize, "%.8g",
1188-
rlinfo->qc_values[shapeindex * rlinfo->band_count + iValue]);
1193+
1194+
double pixelValue =
1195+
rlinfo->qc_values[shapeindex * rlinfo->band_count + iValue];
1196+
1197+
switch (rlinfo->eDataType) {
1198+
case GDT_Byte:
1199+
case GDT_Int8:
1200+
case GDT_UInt16:
1201+
case GDT_Int16:
1202+
case GDT_Int32:
1203+
snprintf(szWork, bufferSize, "%d", (int)pixelValue);
1204+
break;
1205+
case GDT_UInt32:
1206+
snprintf(szWork, bufferSize, "%u", (unsigned int)pixelValue);
1207+
break;
1208+
case GDT_Float64:
1209+
snprintf(szWork, bufferSize, "%.12g", pixelValue);
1210+
break;
1211+
default:
1212+
// GDT_Float32
1213+
snprintf(szWork, bufferSize, "%.8g", pixelValue);
1214+
break;
1215+
}
1216+
11891217
} else if (EQUAL(layer->items[i], "class") && rlinfo->qc_class) {
11901218
int p_class = rlinfo->qc_class[shapeindex];
11911219
if (layer->class[p_class] -> name != NULL)

0 commit comments

Comments
 (0)