/**********************************************************************

 * $Id: mapuv.c 12629 2011-10-06 18:06:34Z aboudreault $

 *

 * Project:  MapServer

 * Purpose:  UV Layer

 * Author:   Alan Boudreault (aboudreault@mapgears.com)

 *

 **********************************************************************

 * Copyright (c) 2011, Alan Boudreault, MapGears

 *

 * Permission is hereby granted, free of charge, to any person obtaining a

 * copy of this software and associated documentation files (the "Software"),

 * to deal in the Software without restriction, including without limitation

 * the rights to use, copy, modify, merge, publish, distribute, sublicense,

 * and/or sell copies of the Software, and to permit persons to whom the

 * Software is furnished to do so, subject to the following conditions:

 *

 * The above copyright notice and this permission notice shall be included in

 * all copies of this Software or works derived from this Software.

 *

 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL

 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING

 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER

 * DEALINGS IN THE SOFTWARE.

 **********************************************************************/

#include "mapserver.h"

#include <assert.h>

#include <math.h>

#include <stdbool.h>

#include "mapows.h"

#include "mapresample.h"

#include "mapthread.h"

#include <cmath>

#include <limits>

#define MSUVRASTER_NUMITEMS 6

#define MSUVRASTER_ANGLE "uv_angle"

#define MSUVRASTER_ANGLEINDEX -100

#define MSUVRASTER_MINUS_ANGLE "uv_minus_angle"

#define MSUVRASTER_MINUSANGLEINDEX -101

#define MSUVRASTER_LENGTH "uv_length"

#define MSUVRASTER_LENGTHINDEX -102

#define MSUVRASTER_LENGTH_2 "uv_length_2"

#define MSUVRASTER_LENGTH2INDEX -103

#define MSUVRASTER_U "u"

#define MSUVRASTER_UINDEX -104

#define MSUVRASTER_V "v"

#define MSUVRASTER_VINDEX -105

#define MSUVRASTER_LON "lon"

#define MSUVRASTER_LONINDEX -106

#define MSUVRASTER_LAT "lat"

#define MSUVRASTER_LATINDEX -107

struct uvRasterLayerInfo {

  /* query cache results */

  int query_results = 0;

  int refcount = 0;

  float *u = nullptr; /* u values */

  float *v = nullptr; /* v values */

  int width = 0;

  int height = 0;

  rectObj extent{};

  int next_shape = 0;

  /* To improve performance of GetShape() when queried on increasing shapeindex

   */

  long last_queried_shapeindex = 0; // value in [0, query_results[ range

  size_t last_raster_off = 0;       // value in [0, width*height[ range

  bool needsLonLat = false;

  reprojectionObj *reprojectorToLonLat = nullptr;

  /* set if the map->extent and map->projection are

     valid in msUVRASTERLayerWhichShapes() */

  mapObj *mapToUseForWhichShapes = nullptr;

  std::string timestring{};

  std::string timefield{};

};

static uvRasterLayerInfo *getLayerInfo(layerObj *layer) {

  return static_cast<uvRasterLayerInfo *>(layer->layerinfo);

}

void msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes(layerObj *layer,

                                                                mapObj *map) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  uvlinfo->mapToUseForWhichShapes = map;

}

static int msUVRASTERLayerInitItemInfo(layerObj *layer) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  int i;

  int *itemindexes;

  int failed = 0;

  if (layer->numitems == 0)

    return MS_SUCCESS;

  if (uvlinfo == NULL) {

    msSetError(MS_MISCERR, "Assertion failed: GDAL layer not opened!!!",

               "msUVRASTERLayerInitItemInfo()");

    return (MS_FAILURE);

  }

  if (layer->iteminfo)

    free(layer->iteminfo);

  if ((layer->iteminfo = (int *)malloc(sizeof(int) * layer->numitems)) ==

      NULL) {

    msSetError(MS_MEMERR, NULL, "msUVRASTERLayerInitItemInfo()");

    return (MS_FAILURE);

  }

  itemindexes = (int *)layer->iteminfo;

  for (i = 0; i < layer->numitems; i++) {

    /* Special attribute names. */

    if (EQUAL(layer->items[i], MSUVRASTER_ANGLE))

      itemindexes[i] = MSUVRASTER_ANGLEINDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_MINUS_ANGLE))

      itemindexes[i] = MSUVRASTER_MINUSANGLEINDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_LENGTH))

      itemindexes[i] = MSUVRASTER_LENGTHINDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_LENGTH_2))

      itemindexes[i] = MSUVRASTER_LENGTH2INDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_U))

      itemindexes[i] = MSUVRASTER_UINDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_V))

      itemindexes[i] = MSUVRASTER_VINDEX;

    else if (EQUAL(layer->items[i], MSUVRASTER_LON)) {

      uvlinfo->needsLonLat = true;

      itemindexes[i] = MSUVRASTER_LONINDEX;

    } else if (EQUAL(layer->items[i], MSUVRASTER_LAT)) {

      uvlinfo->needsLonLat = true;

      itemindexes[i] = MSUVRASTER_LATINDEX;

    } else {

      itemindexes[i] = -1;

      msSetError(MS_MISCERR, "Invalid Field name: %s",

                 "msUVRASTERLayerInitItemInfo()", layer->items[i]);

      failed = 1;

    }

  }

  return failed ? (MS_FAILURE) : (MS_SUCCESS);

}

void msUVRASTERLayerFreeItemInfo(layerObj *layer) {

  if (layer->iteminfo)

    free(layer->iteminfo);

  layer->iteminfo = NULL;

}

static void msUVRasterLayerInfoInitialize(layerObj *layer) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  if (uvlinfo != NULL)

    return;

  uvlinfo = new uvRasterLayerInfo;

  layer->layerinfo = uvlinfo;

  /* Set attribute type to Real, unless the user has explicitly set */

  /* something else. */

  {

    const char *const items[] = {

        MSUVRASTER_ANGLE,    MSUVRASTER_MINUS_ANGLE, MSUVRASTER_LENGTH,

        MSUVRASTER_LENGTH_2, MSUVRASTER_U,           MSUVRASTER_V,

    };

    size_t i;

    for (i = 0; i < sizeof(items) / sizeof(items[0]); ++i) {

      char szTmp[100];

      snprintf(szTmp, sizeof(szTmp), "%s_type", items[i]);

      if (msOWSLookupMetadata(&(layer->metadata), "OFG", szTmp) == NULL) {

        snprintf(szTmp, sizeof(szTmp), "gml_%s_type", items[i]);

        msInsertHashTable(&(layer->metadata), szTmp, "Real");

      }

    }

  }

}

static void msUVRasterLayerInfoFree(layerObj *layer)

{

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  if (uvlinfo == NULL)

    return;

  free(uvlinfo->u);

  free(uvlinfo->v);

  if (uvlinfo->reprojectorToLonLat) {

    msProjectDestroyReprojector(uvlinfo->reprojectorToLonLat);

  }

  delete uvlinfo;

  layer->layerinfo = NULL;

}

int msUVRASTERLayerOpen(layerObj *layer) {

  /* If we don't have info, initialize an empty one now */

  if (layer->layerinfo == NULL)

    msUVRasterLayerInfoInitialize(layer);

  if (layer->layerinfo == NULL)

    return MS_FAILURE;

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  uvlinfo->refcount = uvlinfo->refcount + 1;

  return MS_SUCCESS;

}

int msUVRASTERLayerIsOpen(layerObj *layer) {

  if (layer->layerinfo)

    return MS_TRUE;

  return MS_FALSE;

}

int msUVRASTERLayerClose(layerObj *layer) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  if (uvlinfo != NULL) {

    uvlinfo->refcount--;

    if (uvlinfo->refcount < 1)

      msUVRasterLayerInfoFree(layer);

  }

  return MS_SUCCESS;

}

int msUVRASTERLayerGetItems(layerObj *layer) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  if (uvlinfo == NULL)

    return MS_FAILURE;

  layer->numitems = 0;

  layer->items = (char **)msSmallCalloc(sizeof(char *), 10);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_ANGLE);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_MINUS_ANGLE);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LENGTH);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LENGTH_2);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_U);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_V);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LON);

  layer->items[layer->numitems++] = msStrdup(MSUVRASTER_LAT);

  layer->items[layer->numitems] = NULL;

  return msUVRASTERLayerInitItemInfo(layer);

}

/**********************************************************************

 *                     msUVRASTERGetValues()

 *

 * Special attribute names are used to return some UV params: uv_angle,

 * uv_length, u and v.

 **********************************************************************/

static char **msUVRASTERGetValues(layerObj *layer, float u, float v,

                                  const pointObj *point) {

  char **values;

  int i = 0;

  char tmp[100];

  float size_scale;

  int *itemindexes = (int *)layer->iteminfo;

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  double lon = HUGE_VAL;

  double lat = HUGE_VAL;

  if (layer->numitems == 0)

    return (NULL);

  if (!layer->iteminfo) { /* Should not happen... but just in case! */

    if (msUVRASTERLayerInitItemInfo(layer) != MS_SUCCESS)

      return NULL;

    itemindexes = (int *)layer->iteminfo; /* reassign after malloc */

  }

  if ((values = (char **)malloc(sizeof(char *) * layer->numitems)) == NULL) {

    msSetError(MS_MEMERR, NULL, "msUVRASTERGetValues()");

    return (NULL);

  }

  /* -------------------------------------------------------------------- */

  /*    Determine desired size_scale.  Default to 1 if not otherwise set  */

  /* -------------------------------------------------------------------- */

  size_scale = 1;

  if (CSLFetchNameValue(layer->processing, "UV_SIZE_SCALE") != NULL) {

    size_scale = atof(CSLFetchNameValue(layer->processing, "UV_SIZE_SCALE"));

  }

  if (uvlinfo->needsLonLat) {

    if (uvlinfo->reprojectorToLonLat == NULL)

      uvlinfo->reprojectorToLonLat =

          msProjectCreateReprojector(&layer->projection, NULL);

    if (uvlinfo->reprojectorToLonLat) {

      pointObj pointWrk = *point;

      if (msProjectPointEx(uvlinfo->reprojectorToLonLat, &pointWrk) ==

          MS_SUCCESS) {

        lon = pointWrk.x;

        lat = pointWrk.y;

      }

    }

  }

  for (i = 0; i < layer->numitems; i++) {

    if (itemindexes[i] == MSUVRASTER_ANGLEINDEX) {

      snprintf(tmp, 100, "%f", (atan2((double)v, (double)u) * 180 / MS_PI));

      values[i] = msStrdup(tmp);

    } else if (itemindexes[i] == MSUVRASTER_MINUSANGLEINDEX) {

      double minus_angle;

      minus_angle = (atan2((double)v, (double)u) * 180 / MS_PI) + 180;

      if (minus_angle >= 360)

        minus_angle -= 360;

      snprintf(tmp, 100, "%f", minus_angle);

      values[i] = msStrdup(tmp);

    } else if ((itemindexes[i] == MSUVRASTER_LENGTHINDEX) ||

               (itemindexes[i] == MSUVRASTER_LENGTH2INDEX)) {

      float length = sqrt((u * u) + (v * v)) * size_scale;

      if (itemindexes[i] == MSUVRASTER_LENGTHINDEX)

        snprintf(tmp, 100, "%f", length);

      else

        snprintf(tmp, 100, "%f", length / 2);

      values[i] = msStrdup(tmp);

    } else if (itemindexes[i] == MSUVRASTER_UINDEX) {

      snprintf(tmp, 100, "%f", u);

      values[i] = msStrdup(tmp);

    } else if (itemindexes[i] == MSUVRASTER_VINDEX) {

      snprintf(tmp, 100, "%f", v);

      values[i] = msStrdup(tmp);

    } else if (itemindexes[i] == MSUVRASTER_LONINDEX) {

      snprintf(tmp, sizeof(tmp), "%.18g", lon);

      values[i] = msStrdup(tmp);

    } else if (itemindexes[i] == MSUVRASTER_LATINDEX) {

      snprintf(tmp, sizeof(tmp), "%.18g", lat);

      values[i] = msStrdup(tmp);

    } else {

      values[i] = NULL;

    }

  }

  return values;

}

rectObj msUVRASTERGetSearchRect(layerObj *layer, mapObj *map) {

  rectObj searchrect = map->extent;

  int bDone = MS_FALSE;

  /* For UVRaster, it is important that the searchrect is not too large */

  /* to avoid insufficient intermediate raster resolution, which could */

  /* happen if we use the default code path, given potential reprojection */

  /* issues when using a map extent that is not in the validity area of */

  /* the layer projection. */

  if (!layer->projection.gt.need_geotransform &&

      !(msProjIsGeographicCRS(&(map->projection)) &&

        msProjIsGeographicCRS(&(layer->projection)))) {

    rectObj layer_ori_extent;

    if (msLayerGetExtent(layer, &layer_ori_extent) == MS_SUCCESS) {

      projectionObj map_proj;

      double map_extent_minx = map->extent.minx;

      double map_extent_miny = map->extent.miny;

      double map_extent_maxx = map->extent.maxx;

      double map_extent_maxy = map->extent.maxy;

      rectObj layer_extent = layer_ori_extent;

      /* Create a variant of map->projection without geotransform for */

      /* conveniency */

      msInitProjection(&map_proj);

      msCopyProjection(&map_proj, &map->projection);

      map_proj.gt.need_geotransform = MS_FALSE;

      if (map->projection.gt.need_geotransform) {

        map_extent_minx =

            map->projection.gt.geotransform[0] +

            map->projection.gt.geotransform[1] * map->extent.minx +

            map->projection.gt.geotransform[2] * map->extent.miny;

        map_extent_miny =

            map->projection.gt.geotransform[3] +

            map->projection.gt.geotransform[4] * map->extent.minx +

            map->projection.gt.geotransform[5] * map->extent.miny;

        map_extent_maxx =

            map->projection.gt.geotransform[0] +

            map->projection.gt.geotransform[1] * map->extent.maxx +

            map->projection.gt.geotransform[2] * map->extent.maxy;

        map_extent_maxy =

            map->projection.gt.geotransform[3] +

            map->projection.gt.geotransform[4] * map->extent.maxx +

            map->projection.gt.geotransform[5] * map->extent.maxy;

      }

      /* Reproject layer extent to map projection */

      msProjectRect(&layer->projection, &map_proj, &layer_extent);

      if (layer_extent.minx <= map_extent_minx &&

          layer_extent.miny <= map_extent_miny &&

          layer_extent.maxx >= map_extent_maxx &&

          layer_extent.maxy >= map_extent_maxy) {

        /* do nothing special if area to map is inside layer extent */

      } else {

        if (layer_extent.minx >= map_extent_minx &&

            layer_extent.maxx <= map_extent_maxx &&

            layer_extent.miny >= map_extent_miny &&

            layer_extent.maxy <= map_extent_maxy) {

          /* if the area to map is larger than the layer extent, then */

          /* use full layer extent and add some margin to reflect the */

          /* proportion of the useful area over the requested bbox */

          double extra_x = (map_extent_maxx - map_extent_minx) /

                           (layer_extent.maxx - layer_extent.minx) *

                           (layer_ori_extent.maxx - layer_ori_extent.minx);

          double extra_y = (map_extent_maxy - map_extent_miny) /

                           (layer_extent.maxy - layer_extent.miny) *

                           (layer_ori_extent.maxy - layer_ori_extent.miny);

          searchrect.minx = layer_ori_extent.minx - extra_x / 2;

          searchrect.maxx = layer_ori_extent.maxx + extra_x / 2;

          searchrect.miny = layer_ori_extent.miny - extra_y / 2;

          searchrect.maxy = layer_ori_extent.maxy + extra_y / 2;

        } else {

          /* otherwise clip the map extent with the reprojected layer */

          /* extent */

          searchrect.minx = MS_MAX(map_extent_minx, layer_extent.minx);

          searchrect.maxx = MS_MIN(map_extent_maxx, layer_extent.maxx);

          searchrect.miny = MS_MAX(map_extent_miny, layer_extent.miny);

          searchrect.maxy = MS_MIN(map_extent_maxy, layer_extent.maxy);

          /* and reproject into the layer projection */

          msProjectRect(&map_proj, &layer->projection, &searchrect);

        }

        bDone = MS_TRUE;

      }

      msFreeProjection(&map_proj);

    }

  }

  if (!bDone)

    msProjectRect(&map->projection, &layer->projection,

                  &searchrect); /* project the searchrect to source coords */

  return searchrect;

}

int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery) {

  uvRasterLayerInfo *uvlinfo = getLayerInfo(layer);

  imageObj *image_tmp;

  outputFormatObj *outputformat = NULL;

  mapObj *map_tmp;

  double map_cellsize;

  unsigned int spacing;

  int width, height;

  char **alteredProcessing = NULL, *saved_layer_mask;

  char **savedProcessing = NULL;

  int bHasLonWrap = MS_FALSE;

  double dfLonWrap = 0.0;

  rectObj oldLayerExtent;

  char *oldLayerData = NULL;

  projectionObj oldLayerProjection;

  int ret;

  memset(&oldLayerExtent, 0, sizeof(oldLayerExtent));

  memset(&oldLayerProjection, 0, sizeof(oldLayerProjection));

  if (layer->debug)

    msDebug("Entering msUVRASTERLayerWhichShapes().\n");

  if (uvlinfo == NULL)

    return MS_FAILURE;

  if (CSLFetchNameValue(layer->processing, "BANDS") == NULL) {

    msSetError(MS_MISCERR,

               "BANDS processing option is required for UV layer. You have to "

               "specified 2 bands.",

               "msUVRASTERLayerWhichShapes()");

    return MS_FAILURE;

  }

  /*

  ** Allocate mapObj structure

  */

  map_tmp = (mapObj *)msSmallCalloc(sizeof(mapObj), 1);

  if (initMap(map_tmp) == -1) { /* initialize this map */

    msFree(map_tmp);

    return (MS_FAILURE);

  }

  /* -------------------------------------------------------------------- */

  /*      Determine desired spacing.  Default to 32 if not otherwise set  */

  /* -------------------------------------------------------------------- */

  spacing = 32;

  if (CSLFetchNameValue(layer->processing, "UV_SPACING") != NULL) {

    spacing = atoi(CSLFetchNameValue(layer->processing, "UV_SPACING"));

    if (spacing == 0)

      spacing = 32;

  }

  width = (int)(layer->map->width / spacing);

  height = (int)(layer->map->height / spacing);

  /* Initialize our dummy map */

  MS_INIT_COLOR(map_tmp->imagecolor, 255, 255, 255, 255);

  map_tmp->resolution = layer->map->resolution;

  map_tmp->defresolution = layer->map->defresolution;

  outputformat = (outputFormatObj *)msSmallCalloc(1, sizeof(outputFormatObj));

  outputformat->bands = 2;

  outputformat->name = NULL;

  outputformat->driver = NULL;

  outputformat->refcount = 0;

  outputformat->vtable = NULL;

  outputformat->device = NULL;

  outputformat->renderer = MS_RENDER_WITH_RAWDATA;

  outputformat->imagemode = MS_IMAGEMODE_FLOAT32;

  msAppendOutputFormat(map_tmp, outputformat);

  msCopyHashTable(&map_tmp->configoptions, &layer->map->configoptions);

  map_tmp->mappath = msStrdup(layer->map->mappath);

  map_tmp->shapepath = msStrdup(layer->map->shapepath);

  map_tmp->gt.rotation_angle = 0.0;

  /* Custom msCopyProjection() that removes lon_wrap parameter */

  {

    int i;

    map_tmp->projection.numargs = 0;

    map_tmp->projection.gt = layer->projection.gt;

    for (i = 0; i < layer->projection.numargs; i++) {

      if (strncmp(layer->projection.args[i],

                  "lon_wrap=", strlen("lon_wrap=")) == 0) {

        bHasLonWrap = MS_TRUE;

        dfLonWrap = atof(layer->projection.args[i] + strlen("lon_wrap="));

      } else {

        map_tmp->projection.args[map_tmp->projection.numargs++] =

            msStrdup(layer->projection.args[i]);

      }

    }

    if (map_tmp->projection.numargs != 0) {

      msProcessProjection(&(map_tmp->projection));

    }

    map_tmp->projection.wellknownprojection =

        layer->projection.wellknownprojection;

  }

  /* Very special case to improve quality for rasters referenced from lon=0 to

   * 360 */

  /* We create a temporary VRT that swiches the 2 hemispheres, and then we */

  /* modify the georeferencing to be in the more standard [-180, 180] range */

  /* and we adjust the layer->data, extent and projection accordingly */

  if (layer->tileindex == NULL && uvlinfo->mapToUseForWhichShapes &&

      bHasLonWrap && dfLonWrap == 180.0) {

    rectObj layerExtent;

    msLayerGetExtent(layer, &layerExtent);

    if (layerExtent.minx == 0 && layerExtent.maxx == 360) {

      GDALDatasetH hDS = NULL;

      char *decrypted_path;

      if (strncmp(layer->data, "<VRTDataset", strlen("<VRTDataset")) == 0) {

        decrypted_path = msStrdup(layer->data);

      } else {

        char szPath[MS_MAXPATHLEN];

        msTryBuildPath3(szPath, layer->map->mappath, layer->map->shapepath,

                        layer->data);

        decrypted_path = msDecryptStringTokens(layer->map, szPath);

      }

      if (decrypted_path) {

        char **connectionoptions;

        GDALAllRegister();

        connectionoptions =

            msGetStringListFromHashTable(&(layer->connectionoptions));

        hDS = GDALOpenEx(decrypted_path, GDAL_OF_RASTER, NULL,

                         (const char *const *)connectionoptions, NULL);

        CSLDestroy(connectionoptions);

      }

      if (hDS != NULL) {

        int iBand;

        int nXSize = GDALGetRasterXSize(hDS);

        int nYSize = GDALGetRasterYSize(hDS);

        int nBands = GDALGetRasterCount(hDS);

        int nMaxLen = 100 + nBands * (800 + 2 * strlen(decrypted_path));

        int nOffset = 0;

        char *pszInlineVRT = static_cast<char *>(msSmallMalloc(nMaxLen));

        snprintf(pszInlineVRT, nMaxLen,

                 "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">", nXSize,

                 nYSize);

        nOffset = strlen(pszInlineVRT);

        for (iBand = 1; iBand <= nBands; iBand++) {

          const char *pszDataType = "Byte";

          switch (GDALGetRasterDataType(GDALGetRasterBand(hDS, iBand))) {

          case GDT_Byte:

            pszDataType = "Byte";

            break;

          case GDT_Int16:

            pszDataType = "Int16";

            break;

          case GDT_UInt16:

            pszDataType = "UInt16";

            break;

          case GDT_Int32:

            pszDataType = "Int32";

            break;

          case GDT_UInt32:

            pszDataType = "UInt32";

            break;

          case GDT_Float32:

            pszDataType = "Float32";

            break;

          case GDT_Float64:

            pszDataType = "Float64";

            break;

          default:

            break;

          }

          snprintf(pszInlineVRT + nOffset, nMaxLen - nOffset,

                   "    <VRTRasterBand dataType=\"%s\" band=\"%d\">"

                   "        <SimpleSource>"

                   "            <SourceFilename "

                   "relativeToVrt=\"1\"><![CDATA[%s]]></SourceFilename>"

                   "            <SourceBand>%d</SourceBand>"

                   "            <SrcRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" "

                   "ySize=\"%d\"/>"

                   "            <DstRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" "

                   "ySize=\"%d\"/>"

                   "        </SimpleSource>"

                   "        <SimpleSource>"

                   "            <SourceFilename "

                   "relativeToVrt=\"1\"><![CDATA[%s]]></SourceFilename>"

                   "            <SourceBand>%d</SourceBand>"

                   "            <SrcRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" "

                   "ySize=\"%d\"/>"

                   "            <DstRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" "

                   "ySize=\"%d\"/>"

                   "        </SimpleSource>"

                   "    </VRTRasterBand>",

                   pszDataType, iBand, decrypted_path, iBand, nXSize / 2, 0,

                   nXSize - nXSize / 2, nYSize, 0, 0, nXSize - nXSize / 2,

                   nYSize, decrypted_path, iBand, 0, 0, nXSize / 2, nYSize,

                   nXSize - nXSize / 2, 0, nXSize / 2, nYSize);

          nOffset += strlen(pszInlineVRT + nOffset);

        }

        snprintf(pszInlineVRT + nOffset, nMaxLen - nOffset, "</VRTDataset>");

        oldLayerExtent = layer->extent;

        oldLayerData = layer->data;

        oldLayerProjection = layer->projection;

        layer->extent.minx = -180;

        layer->extent.maxx = 180;

        layer->data = pszInlineVRT;

        layer->projection = map_tmp->projection;

        /* map_tmp->projection is actually layer->projection without lon_wrap */

        rect = uvlinfo->mapToUseForWhichShapes->extent;

        msProjectRect(&uvlinfo->mapToUseForWhichShapes->projection,

                      &map_tmp->projection, &rect);

        bHasLonWrap = MS_FALSE;

        GDALClose(hDS);

      }

      msFree(decrypted_path);

    }

  }

  if (isQuery) {

    /* For query mode, use layer->map->extent reprojected rather than */

    /* the provided rect. Generic query code will filter returned features. */

    rect = msUVRASTERGetSearchRect(layer, layer->map);

  }

  map_cellsize = MS_MAX(MS_CELLSIZE(rect.minx, rect.maxx, layer->map->width),

                        MS_CELLSIZE(rect.miny, rect.maxy, layer->map->height));

  map_tmp->cellsize = map_cellsize * spacing;

  map_tmp->extent.minx =

      rect.minx - (0.5 * map_cellsize) + (0.5 * map_tmp->cellsize);

  map_tmp->extent.miny =

      rect.miny - (0.5 * map_cellsize) + (0.5 * map_tmp->cellsize);

  map_tmp->extent.maxx =

      map_tmp->extent.minx + ((width - 1) * map_tmp->cellsize);

  map_tmp->extent.maxy =

      map_tmp->extent.miny + ((height - 1) * map_tmp->cellsize);

  if (bHasLonWrap && dfLonWrap == 180.0) {

    if (map_tmp->extent.minx >= 180) {

      /* Request on the right half of the shifted raster (= western hemisphere)

       */

      map_tmp->extent.minx -= 360;

      map_tmp->extent.maxx -= 360;

    } else if (map_tmp->extent.maxx >= 180.0) {

      /* Request spanning on the 2 hemispheres => drawing whole planet */

      /* Take only into account vertical resolution, as horizontal one */

      /* will be unreliable (assuming square pixels...) */

      map_cellsize = MS_CELLSIZE(rect.miny, rect.maxy, layer->map->height);

      map_tmp->cellsize = map_cellsize * spacing;

      width = 360.0 / map_tmp->cellsize;

      map_tmp->extent.minx = -180.0 + (0.5 * map_tmp->cellsize);

      map_tmp->extent.maxx = 180.0 - (0.5 * map_tmp->cellsize);

    }

  }

  if (layer->debug)

    msDebug(

        "msUVRASTERLayerWhichShapes(): width: %d, height: %d, cellsize: %g\n",

        width, height, map_tmp->cellsize);

  if (layer->debug == MS_DEBUGLEVEL_VVV)

    msDebug("msUVRASTERLayerWhichShapes(): extent: %g %g %g %g\n",

            map_tmp->extent.minx, map_tmp->extent.miny, map_tmp->extent.maxx,

            map_tmp->extent.maxy);

  /* important to use that function, to compute map

     geotransform, used by the resampling*/

  msMapSetSize(map_tmp, width, height);

  if (layer->debug == MS_DEBUGLEVEL_VVV)

    msDebug("msUVRASTERLayerWhichShapes(): geotransform: %g %g %g %g %g %g\n",

            map_tmp->gt.geotransform[0], map_tmp->gt.geotransform[1],

            map_tmp->gt.geotransform[2], map_tmp->gt.geotransform[3],

            map_tmp->gt.geotransform[4], map_tmp->gt.geotransform[5]);

  uvlinfo->extent = map_tmp->extent;

  image_tmp = msImageCreate(width, height, map_tmp->outputformatlist[0], NULL,

                            NULL, map_tmp->resolution, map_tmp->defresolution,

                            &(map_tmp->imagecolor));

  /* Default set to AVERAGE resampling */

  if (CSLFetchNameValue(layer->processing, "RESAMPLE") == NULL) {

    alteredProcessing = CSLDuplicate(layer->processing);

    alteredProcessing =

        CSLSetNameValue(alteredProcessing, "RESAMPLE", "AVERAGE");

    savedProcessing = layer->processing;

    layer->processing = alteredProcessing;

  }

  /* disable masking at this level: we don't want to apply the mask at the

   * raster level, it will be applied with the correct cellsize and image size

   * in the vector rendering phase.

   */

  saved_layer_mask = layer->mask;

  layer->mask = NULL;

  if (layer->tileindex) {

    expressionObj old_filter;

    if (!uvlinfo->timestring.empty()) {

      msInitExpression(&old_filter);

      msCopyExpression(&old_filter, &layer->filter); /* save existing filter */

      msFreeExpression(&layer->filter);

      msLayerMakeBackticsTimeFilter(layer, uvlinfo->timestring.c_str(),

                                    uvlinfo->timefield.c_str());

    }

    ret = msDrawRasterLayerLow(map_tmp, layer, image_tmp, NULL);

    if (!uvlinfo->timestring.empty()) {

      msCopyExpression(&layer->filter, &old_filter); /* restore old filter */

      msFreeExpression(&old_filter);

    }

  } else {

    ret = msDrawRasterLayerLow(map_tmp, layer, image_tmp, NULL);

  }

  /* restore layer attributes if we went through the above on-the-fly VRT */

  if (oldLayerData) {

    msFree(layer->data);

    layer->data = oldLayerData;

    layer->extent = oldLayerExtent;

    layer->projection = oldLayerProjection;

  }

  /* restore layer mask */

  layer->mask = saved_layer_mask;

  /* restore the saved processing */

  if (alteredProcessing != NULL) {

    layer->processing = savedProcessing;

    CSLDestroy(alteredProcessing);

  }

  if (ret == MS_FAILURE) {

    msSetError(MS_MISCERR, "Unable to draw raster data.",

               "msUVRASTERLayerWhichShapes()");

    msFreeMap(map_tmp);

    msFreeImage(image_tmp);

    return MS_FAILURE;

  }

  /* free old query arrays */

  free(uvlinfo->u);

  free(uvlinfo->v);

  /* Update our uv layer structure */

  uvlinfo->width = width;

  uvlinfo->height = height;

  uvlinfo->query_results = 0;

  uvlinfo->last_queried_shapeindex = 0;

  uvlinfo->last_raster_off = 0;

  uvlinfo->u = (float *)msSmallMalloc(sizeof(float *) * width * height);

  uvlinfo->v = (float *)msSmallMalloc(sizeof(float *) * width * height);

  reprojectionObj *reprojectorLayerToMap = nullptr;

  if (layer->project) {

    reprojectorLayerToMap =

        msProjectCreateReprojector(&layer->projection, &layer->map->projection);

  }

  // To reproject the wind direction, use a ~ 10 m long line. This is a value

  // not too small and not too big to be able to reproject a small vector.

  double typicalLength = 10.0;

  if (msProjIsGeographicCRS(&(layer->projection))) {

    const double dfDegToMeter =

        msProjGetSemiMajorAxis(&(layer->projection)) * (MS_PI / 180);

    const double dfMeterToDeg = 1.0 / dfDegToMeter;

    typicalLength *= dfMeterToDeg;

  }

  for (size_t off = 0; off < static_cast<size_t>(width) * height; ++off) {

    /* Ignore invalid pixels (at nodata), or (u,v)=(0,0) */

    if (MS_GET_BIT(image_tmp->img_mask, off)) {

      uvlinfo->u[off] = image_tmp->img.raw_float[off];

      uvlinfo->v[off] =

          image_tmp->img.raw_float[off + static_cast<size_t>(width) * height];

      if (!(uvlinfo->u[off] == 0 && uvlinfo->v[off] == 0)) {

        uvlinfo->query_results++;

        if (reprojectorLayerToMap) {

          // Compute coordinates of the pixel in layer projection

          const int x = static_cast<int>(off % uvlinfo->width);

          const int y = static_cast<int>(off / uvlinfo->width);

          const double xLayer =

              Pix2Georef(x, 0, uvlinfo->width - 1, uvlinfo->extent.minx,

                         uvlinfo->extent.maxx, MS_FALSE);

          const double yLayer =

              Pix2Georef(y, 0, uvlinfo->height - 1, uvlinfo->extent.miny,

                         uvlinfo->extent.maxy, MS_TRUE);

          // Creates a vector starting at point.x,point.y, of length

          // ~ 10 m, using the angle defined by (u,v)

          pointObj point1, point2;

          point1.x = xLayer;

          point1.y = yLayer;

          const double lengthInLayerProj =

              hypotf(uvlinfo->u[off], uvlinfo->v[off]);

          point2.x =

              xLayer + typicalLength * (uvlinfo->u[off] / lengthInLayerProj);

          point2.y =

              yLayer + typicalLength * (uvlinfo->v[off] / lengthInLayerProj);

          // Reprojects that vector to map projection

          if (msProjectPointEx(reprojectorLayerToMap, &point1) == MS_SUCCESS &&

              msProjectPointEx(reprojectorLayerToMap, &point2) == MS_SUCCESS) {

            // Now computes the (u,v) component in map projection,

            // preserving its original length

            const double dx = point2.x - point1.x;

            const double dy = point2.y - point1.y;

            const double lengthInMapProj = hypot(dx, dy);