Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
UVRaster: better result quality for rasters whose longitude span from… 
… 0 to 360 deg (complementary fix for #5502)
  • Loading branch information
@rouault
rouault committed Nov 17, 2017
1 parent 4c5fa18 commit 116c7d2
Show file tree
Hide file tree
Showing 8 changed files with 177 additions and 21 deletions.
14 changes: 14 additions & 0 deletions mapdraw.c
View file
Expand Up @@ -959,6 +959,14 @@ int msDrawVectorLayer(mapObj *map, layerObj *layer, imageObj *image)
if((map->projection.numargs > 0) && (layer->projection.numargs > 0)) {
int bDone = MS_FALSE;

if( layer->connectiontype == MS_UVRASTER )
{
/* Nasty hack to make msUVRASTERLayerWhichShapes() aware that the */
/* original area of interest is (map->extent, map->projection)... */
/* Useful when dealin with UVRASTER that extend beyond 180 deg */
msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes( layer, map );
}

/* 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 */
Expand Down Expand Up @@ -1058,6 +1066,12 @@ int msDrawVectorLayer(mapObj *map, layerObj *layer, imageObj *image)
}

status = msLayerWhichShapes(layer, searchrect, MS_FALSE);

if( layer->connectiontype == MS_UVRASTER )
{
msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes( layer, NULL );
}

if(status == MS_DONE) { /* no overlap */
msLayerClose(layer);
return MS_SUCCESS;
Expand Down
15 changes: 12 additions & 3 deletions mapraster.c
View file
Expand Up @@ -690,14 +690,23 @@ int msDrawRasterLayerLow(mapObj *map, layerObj *layer, imageObj *image,
}

if(layer->connectiontype != MS_KERNELDENSITY) {
char* pszPath;
if(strlen(filename) == 0) continue;

if(layer->debug == MS_TRUE)
msDebug( "msDrawRasterLayerLow(%s): Filename is: %s\n", layer->name, filename);

msDrawRasterBuildRasterPath(map, layer, filename, szPath);
if( strncmp(filename, "<VRTDataset", strlen("<VRTDataset")) == 0 )
{
pszPath = filename;
}
else
{
msDrawRasterBuildRasterPath(map, layer, filename, szPath);
pszPath = szPath;
}
if(layer->debug == MS_TRUE)
msDebug("msDrawRasterLayerLow(%s): Path is: %s\n", layer->name, szPath);
msDebug("msDrawRasterLayerLow(%s): Path is: %s\n", layer->name, pszPath);

/*
** Note: because we do decryption after the above path expansion
Expand All @@ -706,7 +715,7 @@ int msDrawRasterLayerLow(mapObj *map, layerObj *layer, imageObj *image,
** components. But that is mostly ok, since stuff like sde,postgres and
** oracle georaster do not use real paths.
*/
decrypted_path = msDecryptStringTokens( map, szPath );
decrypted_path = msDecryptStringTokens( map, pszPath );
if( decrypted_path == NULL )
return MS_FAILURE;

Expand Down
2 changes: 2 additions & 0 deletions mapserver.h
View file
Expand Up @@ -2503,6 +2503,8 @@ void msPopulateTextSymbolForLabelAndString(textSymbolObj *ts, labelObj *l, char
MS_DLL_EXPORT int msUnionLayerInitializeVirtualTable(layerObj *layer);
MS_DLL_EXPORT void msPluginFreeVirtualTableFactory(void);

void msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes(layerObj* layer, mapObj* map);

/* ==================================================================== */
/* Prototypes for functions in mapdraw.c */
/* ==================================================================== */
Expand Down
167 changes: 149 additions & 18 deletions mapuvraster.c
View file
Expand Up @@ -84,8 +84,17 @@ typedef struct {
int next_shape;
int x, y; /* used internally in msUVRasterLayerNextShape() */

mapObj* mapToUseForWhichShapes; /* set if the map->extent and map->projection are valid in msUVRASTERLayerWhichShapes() */

} uvRasterLayerInfo;

void msUVRASTERLayerUseMapExtentAndProjectionForNextWhichShapes(layerObj* layer,
mapObj* map)
{
uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
uvlinfo->mapToUseForWhichShapes = map;
}

static int msUVRASTERLayerInitItemInfo(layerObj *layer)
{
uvRasterLayerInfo *uvlinfo = (uvRasterLayerInfo *) layer->layerinfo;
Expand Down Expand Up @@ -347,7 +356,11 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
char **savedProcessing = NULL;
int bHasLonWrap = MS_FALSE;
double dfLonWrap = 0.0;

rectObj oldLayerExtent;
char* oldLayerData = NULL;
projectionObj oldLayerProjection;
int ret;

if (layer->debug)
msDebug("Entering msUVRASTERLayerWhichShapes().\n");

Expand Down Expand Up @@ -387,9 +400,6 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)

width = (int)ceil(layer->map->width/spacing);
height = (int)ceil(layer->map->height/spacing);
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;

/* Initialize our dummy map */
MS_INIT_COLOR(map_tmp->imagecolor, 255,255,255,255);
Expand All @@ -410,10 +420,6 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
msCopyHashTable(&map_tmp->configoptions, &layer->map->configoptions);
map_tmp->mappath = msStrdup(layer->map->mappath);
map_tmp->shapepath = msStrdup(layer->map->shapepath);
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);
map_tmp->gt.rotation_angle = 0.0;

/* Custom msCopyProjection() that removes lon_wrap parameter */
Expand Down Expand Up @@ -443,6 +449,124 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
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 georeferncing 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 )
{
GDALAllRegister();
hDS = GDALOpen(decrypted_path, GA_ReadOnly );
}
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 = 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 );
}
}


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) */
Expand Down Expand Up @@ -503,17 +627,15 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
*/
saved_layer_mask = layer->mask;
layer->mask = NULL;
if (msDrawRasterLayerLow(map_tmp, layer, image_tmp, NULL ) == MS_FAILURE) {
msSetError(MS_MISCERR, "Unable to draw raster data.", "msUVRASTERLayerWhichShapes()");
layer->mask = saved_layer_mask;
ret = msDrawRasterLayerLow(map_tmp, layer, image_tmp, NULL );

if (alteredProcessing != NULL) {
layer->processing = savedProcessing;
CSLDestroy(alteredProcessing);
}
msFreeMap(map_tmp);
msFreeImage(image_tmp);
return MS_FAILURE;
/* 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 */
Expand All @@ -525,6 +647,15 @@ int msUVRASTERLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
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 */
if (uvlinfo->u) {
for (i=0; i<uvlinfo->width; ++i) {
Expand Down
Binary file modified msautotest/wxs/expected/wms_uvraster_lonwrap180_accross_hemisphere.png
View file
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified msautotest/wxs/expected/wms_uvraster_lonwrap180_east_hemisphere.png
View file
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified msautotest/wxs/expected/wms_uvraster_lonwrap180_west_hemisphere.png
View file
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified msautotest/wxs/expected/wms_uvraster_lonwrap180_whole_world.png
View file
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit 116c7d2

Please sign in to comment.