Ignore:
Timestamp:
Jun 2, 2017, 6:03:41 PM (7 years ago)
Author:
cresson
Message:

ENH: Enhance raster styling in mapserver support

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/ms-style/zoo-project/zoo-kernel/service_internal_ms.c

    r781 r830  
    3434#include "server_internal.h"
    3535#include "response_print.h"
     36
     37static double PRECISION = 0.001;
     38static int MAX_NUMBER_STRING_SIZE = 32;
     39
     40/**
     41 * Double to ASCII
     42 */
     43char * dtoa(char *s, double n) {
     44
     45  // handle special cases
     46  if (isnan(n)) {
     47    strcpy(s, "nan");
     48  } else if (isinf(n)) {
     49    strcpy(s, "inf");
     50  } else if (n == 0.0) {
     51    strcpy(s, "0");
     52  } else {
     53    int digit, m, m1;
     54    char *c = s;
     55    int neg = (n < 0);
     56    if (neg)
     57      n = -n;
     58    // calculate magnitude
     59    m = log10(n);
     60    int useExp = (m >= 14 || (neg && m >= 9) || m <= -9);
     61    if (neg)
     62      *(c++) = '-';
     63    // set up for scientific notation
     64    if (useExp) {
     65      if (m < 0)
     66        m -= 1.0;
     67      n = n / pow(10.0, m);
     68      m1 = m;
     69      m = 0;
     70    }
     71    if (m < 1.0) {
     72      m = 0;
     73    }
     74    // convert the number
     75    while (n > PRECISION || m >= 0) {
     76      double weight = pow(10.0, m);
     77      if (weight > 0 && !isinf(weight)) {
     78        digit = floor(n / weight);
     79        n -= (digit * weight);
     80        *(c++) = '0' + digit;
     81      }
     82      if (m == 0 && n > 0)
     83        *(c++) = '.';
     84      m--;
     85    }
     86    if (useExp) {
     87      // convert the exponent
     88      int i, j;
     89      *(c++) = 'e';
     90      if (m1 > 0) {
     91        *(c++) = '+';
     92      } else {
     93        *(c++) = '-';
     94        m1 = -m1;
     95      }
     96      m = 0;
     97      while (m1 > 0) {
     98        *(c++) = '0' + m1 % 10;
     99        m1 /= 10;
     100        m++;
     101      }
     102      c -= m;
     103      for (i = 0, j = m-1; i<j; i++, j--) {
     104        // swap without temporary
     105        c[i] ^= c[j];
     106        c[j] ^= c[i];
     107        c[i] ^= c[j];
     108      }
     109      c += m;
     110    }
     111    *(c) = '\0';
     112  }
     113  return s;
     114}
     115
     116/**
     117 * List of allowed raster styles
     118 */
     119enum MS_RASTER_STYLES{
     120  LINEAR_STRETCHING,
     121  COLOR_PALETTE
     122};
     123enum LINEAR_STRETCHING_TYPE{
     124  MINMAX,
     125  MEANSTD
     126};
    36127
    37128/**
    … …  
    727818 */
    728819int tryGdal(maps* conf,maps* output,mapObj* m){
     820
     821  /*
     822   * Detect the raster style
     823   */
     824
     825  /* msRasterResample (NEAREST/AVERAGE/BILINEAR) */
     826  const char * msRasterResamplingType               = "msRasterResample";
     827  /* msRasterStyle (linearStretching/colorPalette) */
     828  const char * msRasterStylePropertyName            = "msRasterStyle";
     829  const char * msRasterStyleLinearStretchingPropertyValue       = "linearStretching";
     830  const char * msRasterStyleColorPalettePropertyValue           = "colorPalette";
     831  /* msRasterStyleOptions (minMax/meanstd) */
     832  const char * msRasterStyleOptionsPropertyName     = "msRasterStyleOptions";
     833  const char * msRasterStyleLinearStretchingMinMaxPropertyName  = "minMax";
     834  const char * msRasterStyleLinearStretchingMeanStdPropertyName = "meanStd";
     835
     836  // Default raster style
     837  int defaultStyleType = LINEAR_STRETCHING;
     838  int defaultLinearstretchingType = MEANSTD;
     839
     840  // Check if there is a defined raster style type
     841  int styleType = defaultStyleType;
     842  int linearStretchingType = defaultLinearstretchingType;
     843  map* msRasterStyle=getMap(output->content, msRasterStylePropertyName);
     844  char * msRasterStyleOptionsContent = "";
     845  if(msRasterStyle!=NULL)
     846    {
     847    // Check if there is options attached
     848    map* msRasterStyleOptions=getMap(output->content, msRasterStyleOptionsPropertyName);
     849    if (msRasterStyleOptions!=NULL)
     850      {
     851      msRasterStyleOptionsContent = msRasterStyleOptions->value;
     852      }
     853
     854    // Detect the raster style
     855#ifdef DEBUGMS
     856    fprintf(stderr,"Detect the raster style %s\n", msRasterStyle->value);
     857#endif
     858    if (strncasecmp(msRasterStyle->value, msRasterStyleLinearStretchingPropertyValue,
     859        strlen(msRasterStyleLinearStretchingPropertyValue))==0)
     860      {
     861#ifdef DEBUGMS
     862      fprintf(stderr,"The raster style is linear stretching\n");
     863#endif
     864      styleType = LINEAR_STRETCHING;
     865      if (strlen(msRasterStyleOptionsContent)>0)
     866        {
     867        if (strncasecmp(msRasterStyleOptionsContent, msRasterStyleLinearStretchingMinMaxPropertyName,
     868            strlen(msRasterStyleLinearStretchingMinMaxPropertyName))==0)
     869          {
     870          linearStretchingType = MINMAX;
     871#ifdef DEBUGMS
     872          fprintf(stderr,"The raster style linear stretching is minmax\n");
     873#endif
     874          }
     875        else if (strncasecmp(msRasterStyleOptionsContent, msRasterStyleLinearStretchingMeanStdPropertyName,
     876            strlen(msRasterStyleLinearStretchingMeanStdPropertyName))==0)
     877          {
     878          linearStretchingType = MEANSTD;
     879#ifdef DEBUGMS
     880          fprintf(stderr,"The raster style linear stretching is meanstd\n");
     881#endif
     882          }
     883        else
     884          {
     885          fprintf(stderr,"Unknown raster style linear stretching method: %s\n", msRasterStyleOptionsContent);
     886          }
     887        }
     888      else
     889        {
     890        fprintf(stderr,"Using default linear stretching type.");
     891        }
     892      }
     893    else if (strncasecmp(msRasterStyle->value, msRasterStyleColorPalettePropertyValue,
     894        strlen(msRasterStyleColorPalettePropertyValue))==0)
     895      {
     896#ifdef DEBUGMS
     897      fprintf(stderr,"The raster style is color palette\n");
     898#endif
     899      styleType = COLOR_PALETTE;
     900      }
     901    else
     902      {
     903      fprintf(stderr,"Unknown raster style: %s. Using default.", msRasterStyle->value);
     904      }
     905    }
     906
     907#ifdef DEBUGMS
     908  fprintf(stderr,"RasterStyle=%i Options=%s\n",styleType,msRasterStyleOptionsContent);
     909#endif
     910
     911  /*
     912   * Get storage
     913   */
    729914  map* tmpMap=getMap(output->content,"storage");
    730915  char *pszFilename=tmpMap->value;
    … …  
    733918  double adfGeoTransform[6];
    734919  int i, iBand;
    735  
     920
    736921  /**
    737922   * Try to open the DataSource using GDAL
    … …  
    791976  m->width=GDALGetRasterXSize( hDataset );
    792977  m->height=GDALGetRasterYSize( hDataset );
    793  
     978
    794979  /**
    795980   * Set projection using Authority Code and Name if available or fallback to
    … …  
    808993  }
    809994
    810 
    811995  /**
    812996   * Set extent
    … …  
    8161000
    8171001      double minX = adfGeoTransform[0]
    818         + adfGeoTransform[2] * GDALGetRasterYSize(hDataset);
     1002                                    + adfGeoTransform[2] * GDALGetRasterYSize(hDataset);
    8191003      double minY = adfGeoTransform[3]
    820         + adfGeoTransform[5] * GDALGetRasterYSize(hDataset);
     1004                                    + adfGeoTransform[5] * GDALGetRasterYSize(hDataset);
    8211005
    8221006      double maxX = adfGeoTransform[0]
    823         + adfGeoTransform[1] * GDALGetRasterXSize(hDataset);
     1007                                    + adfGeoTransform[1] * GDALGetRasterXSize(hDataset);
    8241008      double maxY = adfGeoTransform[3]
    825         + adfGeoTransform[4] * GDALGetRasterXSize(hDataset);
    826 
    827        setMsExtent(output,m,myLayer,minX,minY,maxX,maxY);
     1009                                    + adfGeoTransform[4] * GDALGetRasterXSize(hDataset);
     1010
     1011      setMsExtent(output,m,myLayer,minX,minY,maxX,maxY);
    8281012
    8291013    }
    … …  
    8571041    }else{
    8581042      if(iBand<4){
    859         char *tmpS=zStrdup(nameBands);
    860         nameBands=(char*)realloc(nameBands,(strlen(nameBands)+strlen(lBands)+1)*sizeof(char));
    861         sprintf(nameBands,"%s %s",tmpS,lBands);
    862         free(tmpS);
     1043        char *tmpS=zStrdup(nameBands);
     1044        nameBands=(char*)realloc(nameBands,(strlen(nameBands)+strlen(lBands)+1)*sizeof(char));
     1045        sprintf(nameBands,"%s %s",tmpS,lBands);
     1046        free(tmpS);
    8631047      }
    8641048    }
    8651049  }
    8661050  msInsertHashTable(&(myLayer->metadata), "ows_bandnames", nameBands);
    867  
     1051
    8681052  /**
    8691053   * Loops over metadata information to setup specific information
    8701054   */
    871   for( iBand = 0; iBand < nBandsI; iBand++ ){
    872     //int         bGotNodata;//, bSuccess;
    873     double      adfCMinMax[2]/*, dfNoData*/;
    874     //int         nBlockXSize, nBlockYSize, nMaskFlags;
    875     //double      /*dfMean, dfStdDev*/;
     1055  for( iBand = 0; iBand < nBandsI; iBand++ )
     1056    {
     1057
     1058    // Compute statistics of the current band
    8761059    hBand = GDALGetRasterBand( hDataset, iBand+1 );
    877 
    8781060    CPLErrorReset();
    879     GDALComputeRasterMinMax( hBand, FALSE, adfCMinMax );
    880     char tmpN[21];
    881     sprintf(tmpN,"Band%d",iBand+1);
    882     if (CPLGetLastErrorType() == CE_None){
    883       char tmpMm[100];
    884       sprintf(tmpMm,"%.3f %.3f",adfCMinMax[0],adfCMinMax[1]);
    885       char tmpI[21];
    886       sprintf(tmpI,"%s_interval",tmpN);
    887       msInsertHashTable(&(myLayer->metadata), tmpI, tmpMm);
    888 
    889       map* test=getMap(output->content,"msClassify");
    890       if(test!=NULL && strncasecmp(test->value,"true",4)==0){
    891         /**
    892          * Classify one band raster pixel value using regular interval
    893          */
    894         int _tmpColors[10][3]={
    895           {102,153,204},
    896           {51,102,153},
    897           {102,102,204},
    898           {51,204,0},
    899           {153,255,102},
    900           {204,255,102},
    901           {102,204,153},
    902           {255,69,64},
    903           {255,192,115},
    904           {255,201,115}
    905         };
    906          
    907         if(nBandsI==1){
    908           double delta=adfCMinMax[1]-adfCMinMax[0];
    909           double interval=delta/10;
    910           double cstep=adfCMinMax[0];
    911           for(i=0;i<10;i++){
    912             /**
    913              * Create a new class
    914              */
    915             if(msGrowLayerClasses(myLayer) == NULL)
    916               return -1;
    917             if(initClass((myLayer->CLASS[myLayer->numclasses])) == -1)
    918               return -1;
    919             if(msGrowClassStyles(myLayer->CLASS[myLayer->numclasses]) == NULL)
    920               return -1;
    921             if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
    922               return -1;
    923            
    924             /**
    925              * Set class name
    926              */
    927             char className[7];
    928             sprintf(className,"class%d",i);
    929             myLayer->CLASS[myLayer->numclasses]->name=zStrdup(className);
    930            
    931             /**
    932              * Set expression
    933              */
    934             char expression[1024];
    935             if(i+1<10)
    936               sprintf(expression,"([pixel]>=%.3f AND [pixel]<%.3f)",cstep,cstep+interval);
    937             else
    938               sprintf(expression,"([pixel]>=%.3f AND [pixel]<=%.3f)",cstep,cstep+interval);
    939             msLoadExpressionString(&myLayer->CLASS[myLayer->numclasses]->expression,expression);
    940            
    941             /**
    942              * Set color
    943              */
    944             myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.red=_tmpColors[i][0];
    945             myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.green=_tmpColors[i][1];
    946             myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.blue=_tmpColors[i][2];
    947             cstep+=interval;
    948             myLayer->CLASS[myLayer->numclasses]->numstyles++;
    949             myLayer->numclasses++;
    950            
    951           }
    952          
    953           char tmpMm[100];
    954           sprintf(tmpMm,"%.3f %.3f",adfCMinMax[0],adfCMinMax[1]);
    955          
    956         }
    957       }
    958       else{
    959         if(nBandsI==1){
    960           myLayer->offsite.red=0;
    961           myLayer->offsite.green=0;
    962           myLayer->offsite.blue=0;
    963         }
    964         msLayerAddProcessing(myLayer,"RESAMPLE=BILINEAR");
    965       }
    966     }
    967     if( strlen(GDALGetRasterUnitType(hBand)) > 0 ){
     1061    double min = 0.0;
     1062    double max = 0.0;
     1063    double mean = 0.0;
     1064    double std = 0.0;
     1065    GDALComputeRasterStatistics  (hBand, 1, &min, &max, &mean, &std, NULL, NULL);
     1066
     1067    char bandIdentifier[21];
     1068    sprintf(bandIdentifier,"Band%d",iBand+1);
     1069    if (CPLGetLastErrorType() == CE_None)
     1070      {
     1071      char bandInterval[100];
     1072      sprintf(bandInterval,"%.3f %.3f",min,max);
     1073      char bandIntervalIdentifier[21];
     1074      sprintf(bandIntervalIdentifier,"%s_interval",bandIdentifier);
     1075      msInsertHashTable(&(myLayer->metadata), bandIntervalIdentifier, bandInterval);
     1076
     1077      // Apply the raster style
     1078      if(styleType == LINEAR_STRETCHING)
     1079        {
     1080
     1081        char msProcessingInstruction[1024];
     1082        double low = mean - 2*std;
     1083        double hi = mean + 2*std;
     1084
     1085        char s1[MAX_NUMBER_STRING_SIZE];
     1086        char s2[MAX_NUMBER_STRING_SIZE];
     1087
     1088        int bn = iBand+1;
     1089        if(linearStretchingType==MINMAX)
     1090          {
     1091          sprintf(msProcessingInstruction, "SCALE_%d=%s,%s",
     1092              bn,
     1093              dtoa(s1,min),
     1094              dtoa(s2,max));
     1095          }
     1096        else if (linearStretchingType==MEANSTD)
     1097          {
     1098          sprintf(msProcessingInstruction, "SCALE_%d=%s,%s",
     1099              bn,
     1100              dtoa(s1,low),
     1101              dtoa(s2,hi));
     1102          }
     1103
     1104        msLayerAddProcessing(myLayer,msProcessingInstruction);
     1105
     1106        } // styleType is LINEAR_STRETCHING
     1107      else if( styleType == COLOR_PALETTE )
     1108        {
     1109        /**
     1110         * Classify one band raster pixel value using regular interval
     1111         * TODO: parse msRasterStyleOptionsContent to retrieve the colormap
     1112         */
     1113        int discreteColorMap[10][3]={
     1114            {102,153,204},
     1115            {51,102,153},
     1116            {102,102,204},
     1117            {51,204,0},
     1118            {153,255,102},
     1119            {204,255,102},
     1120            {102,204,153},
     1121            {255,69,64},
     1122            {255,192,115},
     1123            {255,201,115}
     1124        };
     1125
     1126        if(nBandsI==1)
     1127          {
     1128          double delta = max - min;
     1129          double interval = delta / 10;
     1130          double cstep = min;
     1131          for(i=0; i<10; i++)
     1132            {
     1133            /**
     1134             * Create a new class
     1135             */
     1136            if(msGrowLayerClasses(myLayer) == NULL)
     1137              return -1;
     1138            if(initClass((myLayer->CLASS[myLayer->numclasses])) == -1)
     1139              return -1;
     1140            if(msGrowClassStyles(myLayer->CLASS[myLayer->numclasses]) == NULL)
     1141              return -1;
     1142            if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
     1143              return -1;
     1144
     1145            /**
     1146             * Set class name
     1147             */
     1148            char className[7];
     1149            sprintf(className,"class%d",i);
     1150            myLayer->CLASS[myLayer->numclasses]->name=zStrdup(className);
     1151
     1152            /**
     1153             * Set expression
     1154             */
     1155            char expression[1024];
     1156            if(i+1<10)
     1157              sprintf(expression,"([pixel]>=%.3f AND [pixel]<%.3f)",cstep,cstep+interval);
     1158            else
     1159              sprintf(expression,"([pixel]>=%.3f AND [pixel]<=%.3f)",cstep,cstep+interval);
     1160            msLoadExpressionString(&myLayer->CLASS[myLayer->numclasses]->expression,expression);
     1161
     1162            /**
     1163             * Set color
     1164             */
     1165            myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.red=discreteColorMap[i][0];
     1166            myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.green=discreteColorMap[i][1];
     1167            myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.blue=discreteColorMap[i][2];
     1168            cstep+=interval;
     1169            myLayer->CLASS[myLayer->numclasses]->numstyles++;
     1170            myLayer->numclasses++;
     1171
     1172            }
     1173
     1174          char tmpMm[100];
     1175          sprintf(tmpMm,"%.3f %.3f",min,max);
     1176
     1177          }
     1178        } // styleType is COLOR_PALETTE
     1179
     1180      } // If no error with GDAL functions
     1181    else
     1182      {
     1183      fprintf(stderr,"Unable to compute raster statistics!\n");
     1184      }
     1185
     1186    /*
     1187     * Set offsite
     1188     */
     1189    int offsiteR = 0;
     1190    int offsiteG = 0;
     1191    int offsiteB = 0;
     1192    int hasNoData = 0;
     1193    double noDataValue = GDALGetRasterNoDataValue(hBand, &hasNoData);
     1194    if (hasNoData)
     1195      {
     1196      offsiteR = (int) noDataValue;
     1197      offsiteG = (int) noDataValue;
     1198      offsiteB = (int) noDataValue;
     1199      }
     1200    myLayer->offsite.red    = offsiteR;
     1201    myLayer->offsite.green  = offsiteG;
     1202    myLayer->offsite.blue   = offsiteB;
     1203
     1204    /*
     1205     * Insert units
     1206     */
     1207    if( strlen(GDALGetRasterUnitType(hBand)) > 0 )
     1208      {
    9681209      char tmpU[21];
    969       sprintf(tmpU,"%s_band_uom",tmpN);
     1210      sprintf(tmpU,"%s_band_uom",bandIdentifier);
    9701211      msInsertHashTable(&(myLayer->metadata), tmpU, GDALGetRasterUnitType(hBand));
    971     }
    972 
    973   }
     1212      }
     1213
     1214    } // next band
     1215
     1216  /*
     1217   * Check if there is resample option
     1218   */
     1219  char msResampleOptionInstruction[1024];
     1220  char * msRasterResampleOptionContent = "BILINEAR";
     1221  map* msRasterResamplingOption=getMap(output->content, msRasterResamplingType);
     1222  if (msRasterResamplingOption!=NULL)
     1223    {
     1224    msRasterResampleOptionContent = msRasterResamplingOption->value;
     1225    }
     1226  sprintf(msResampleOptionInstruction, "RESAMPLE=%s",msRasterResampleOptionContent);
     1227  msLayerAddProcessing(myLayer, msResampleOptionInstruction);
    9741228
    9751229  m->layerorder[m->numlayers] = m->numlayers;
Note: See TracChangeset for help on using the changeset viewer.

Download now !

Search

Context Navigation

Main Navigation

ZOO Sponsors

http://www.zoo-project.org/trac/chrome/site/img/geolabs-logo.pnghttp://www.zoo-project.org/trac/chrome/site/img/neogeo-logo.png http://www.zoo-project.org/trac/chrome/site/img/apptech-logo.png http://www.zoo-project.org/trac/chrome/site/img/3liz-logo.png http://www.zoo-project.org/trac/chrome/site/img/gateway-logo.png

Become a sponsor !

Knowledge partners

http://www.zoo-project.org/trac/chrome/site/img/ocu-logo.png http://www.zoo-project.org/trac/chrome/site/img/gucas-logo.png http://www.zoo-project.org/trac/chrome/site/img/polimi-logo.png http://www.zoo-project.org/trac/chrome/site/img/fem-logo.png http://www.zoo-project.org/trac/chrome/site/img/supsi-logo.png http://www.zoo-project.org/trac/chrome/site/img/cumtb-logo.png

Become a knowledge partner

Related links

http://zoo-project.org/img/ogclogo.png http://zoo-project.org/img/osgeologo.png