Changeset 832 for branches


Ignore:
Timestamp:
Jun 9, 2017, 12:32:35 PM (8 years ago)
Author:
cresson
Message:

ENH: Add custom mapfile styling option

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r830 r832  
    115115
    116116/**
    117  * List of allowed raster styles
     117 * List of allowed raster styles and sub-styles
    118118 */
    119119enum MS_RASTER_STYLES{
    120120  LINEAR_STRETCHING,
    121   COLOR_PALETTE
     121  CLASSIFY
    122122};
    123123enum LINEAR_STRETCHING_TYPE{
    … …  
    125125  MEANSTD
    126126};
     127enum CLASSES_TYPE{
     128  AUTO,
     129  USER
     130};
     131
     132/*
     133 * Functions producing the 'jet' colormap
     134 */
     135double interpolate( double val, double y0, double x0, double y1, double x1 ) {
     136    return (val-x0)*(y1-y0)/(x1-x0) + y0;
     137}
     138
     139double base( double val ) {
     140    if ( val <= 0.25 ) return 0;
     141    else if ( val <= 0.75 ) return interpolate( val, 0.0, 0.25, 1.0, 0.75 );
     142    else if ( val <= 1.25 ) return 1.0;
     143    else if ( val <= 1.75 ) return interpolate( val, 1.0, 1.25, 0.0, 1.75 );
     144    else return 0.0;
     145}
    127146
    128147/**
    … …  
    764783    if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
    765784      return -1;
     785
    766786    /**
    767787     * Apply msStyle else fallback to the default style
    … …  
    824844
    825845  /* msRasterResample (NEAREST/AVERAGE/BILINEAR) */
    826   const char * msRasterResamplingType               = "msRasterResample";
    827   /* msRasterStyle (linearStretching/colorPalette) */
     846  const char * msRasterResamplingPropertyName       = "msRasterResample";
     847  /* msRasterStyle (linearStretching/classify) */
    828848  const char * msRasterStylePropertyName            = "msRasterStyle";
    829849  const char * msRasterStyleLinearStretchingPropertyValue       = "linearStretching";
    830   const char * msRasterStyleColorPalettePropertyValue           = "colorPalette";
    831   /* msRasterStyleOptions (minMax/meanstd) */
     850  const char * msRasterStyleColorPalettePropertyValue           = "classify";
    832851  const char * msRasterStyleOptionsPropertyName     = "msRasterStyleOptions";
     852  /* options for linear stretching */
    833853  const char * msRasterStyleLinearStretchingMinMaxPropertyName  = "minMax";
    834854  const char * msRasterStyleLinearStretchingMeanStdPropertyName = "meanStd";
     855
     856  const unsigned int msRasterStyleClassifyAutoMaximumNumberOfClasses = 256;
    835857
    836858  // Default raster style
    837859  int defaultStyleType = LINEAR_STRETCHING;
    838860  int defaultLinearstretchingType = MEANSTD;
     861  int defaultClassifyType = AUTO;
    839862
    840863  // Check if there is a defined raster style type
    841864  int styleType = defaultStyleType;
    842865  int linearStretchingType = defaultLinearstretchingType;
     866  int classifyType = defaultClassifyType;
    843867  map* msRasterStyle=getMap(output->content, msRasterStylePropertyName);
    844868  char * msRasterStyleOptionsContent = "";
     869  char * msRasterStyleFileContent = "";
    845870  if(msRasterStyle!=NULL)
    846871    {
     872#ifdef DEBUGMS
     873    fprintf(stderr,"msRasterStyle=%s\n", msRasterStyle->value);
     874#endif
    847875    // Check if there is options attached
    848876    map* msRasterStyleOptions=getMap(output->content, msRasterStyleOptionsPropertyName);
    … …  
    850878      {
    851879      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
     880#ifdef DEBUGMS
     881      fprintf(stderr,"msRasterStyleOptions=%s\n", msRasterStyleOptionsContent);
     882#endif
     883
     884      }
     885
    858886    if (strncasecmp(msRasterStyle->value, msRasterStyleLinearStretchingPropertyValue,
    859887        strlen(msRasterStyleLinearStretchingPropertyValue))==0)
    … …  
    863891#endif
    864892      styleType = LINEAR_STRETCHING;
     893
    865894      if (strlen(msRasterStyleOptionsContent)>0)
    866895        {
    … …  
    870899          linearStretchingType = MINMAX;
    871900#ifdef DEBUGMS
    872           fprintf(stderr,"The raster style linear stretching is minmax\n");
     901          fprintf(stderr,"The raster style linear stretching option is minmax\n");
    873902#endif
    874903          }
    … …  
    878907          linearStretchingType = MEANSTD;
    879908#ifdef DEBUGMS
    880           fprintf(stderr,"The raster style linear stretching is meanstd\n");
     909          fprintf(stderr,"The raster style linear stretching option is meanstd\n");
    881910#endif
    882911          }
    … …  
    885914          fprintf(stderr,"Unknown raster style linear stretching method: %s\n", msRasterStyleOptionsContent);
    886915          }
    887         }
     916        } // raster style options (for linear stretching) are not empty
    888917      else
    889918        {
    890         fprintf(stderr,"Using default linear stretching type.");
     919        fprintf(stderr,"Raster style options for linear stretching are empty. Using default.\n");
    891920        }
    892       }
     921      } // raster style is linear stretching
     922
    893923    else if (strncasecmp(msRasterStyle->value, msRasterStyleColorPalettePropertyValue,
    894924        strlen(msRasterStyleColorPalettePropertyValue))==0)
    895925      {
    896926#ifdef DEBUGMS
    897       fprintf(stderr,"The raster style is color palette\n");
    898 #endif
    899       styleType = COLOR_PALETTE;
    900       }
     927      fprintf(stderr,"The raster style is classify\n");
     928#endif
     929      styleType = CLASSIFY;
     930      if (strlen(msRasterStyleOptionsContent)==0)
     931        {
     932        classifyType = AUTO;
     933#ifdef DEBUGMS
     934        fprintf(stderr,"The raster style classes is set to auto\n");
     935#endif
     936        } // raster style is classify, with automatic styling
     937      else
     938        {
     939        classifyType = USER;
     940#ifdef DEBUGMS
     941        fprintf(stderr,"The raster style classes is user defined: %s\n", msRasterStyleOptionsContent);
     942#endif
     943        }
     944      } // raster style is classes
     945
    901946    else
    902947      {
    903948      fprintf(stderr,"Unknown raster style: %s. Using default.", msRasterStyle->value);
    904949      }
    905     }
    906 
    907 #ifdef DEBUGMS
    908   fprintf(stderr,"RasterStyle=%i Options=%s\n",styleType,msRasterStyleOptionsContent);
     950
     951    } // raster style is not null
     952
     953  /*
     954   * This is just for the backward compatibility and will be deprecated
     955   */
     956  {
     957  map* test=getMap(output->content,"msClassify");
     958  if(test!=NULL && strncasecmp(test->value,"true",4)==0)
     959    {
     960    styleType = CLASSIFY;
     961    classifyType = AUTO;
     962    }
     963  }
     964
     965#ifdef DEBUGMS
     966  fprintf(stderr,"Styling options:\n\tRasterStyle=%i\n\tLinearStretching options=%i\n\tClassify options=%i\n",
     967      styleType,linearStretchingType, classifyType);
    909968#endif
    910969
    … …  
    10641123    double std = 0.0;
    10651124    GDALComputeRasterStatistics  (hBand, 1, &min, &max, &mean, &std, NULL, NULL);
    1066 
     1125#ifdef DEBUGMS
     1126      fprintf(stderr,"Computed raster stats for band %i: min=%.3f max=%.3f mean=%.3f std=%.3f\n",
     1127          iBand, min, max, mean, std);
     1128#endif
    10671129    char bandIdentifier[21];
    10681130    sprintf(bandIdentifier,"Band%d",iBand+1);
    … …  
    10791141        {
    10801142
    1081         char msProcessingInstruction[1024];
    1082         double low = mean - 2*std;
    1083         double hi = mean + 2*std;
     1143        char msProcessingDirective[1024];
     1144        double low = 0.0;
     1145        double hi = 1.0;
    10841146
    10851147        char s1[MAX_NUMBER_STRING_SIZE];
    … …  
    10891151        if(linearStretchingType==MINMAX)
    10901152          {
    1091           sprintf(msProcessingInstruction, "SCALE_%d=%s,%s",
    1092               bn,
    1093               dtoa(s1,min),
    1094               dtoa(s2,max));
     1153          low = min;
     1154          hi = max;
    10951155          }
    10961156        else if (linearStretchingType==MEANSTD)
    10971157          {
    1098           sprintf(msProcessingInstruction, "SCALE_%d=%s,%s",
    1099               bn,
    1100               dtoa(s1,low),
    1101               dtoa(s2,hi));
     1158          low = mean - 2*std;
     1159          hi = mean + 2*std;
    11021160          }
    1103 
    1104         msLayerAddProcessing(myLayer,msProcessingInstruction);
     1161#ifdef DEBUGMS
     1162      fprintf(stderr,"Processing the raster using a stretch btw %.3f and %.3f\n", low, hi);
     1163#endif
     1164        sprintf(msProcessingDirective, "SCALE_%d=%s,%s", bn, dtoa(s1,low), dtoa(s2,hi));
     1165        msLayerAddProcessing(myLayer,msProcessingDirective);
    11051166
    11061167        } // styleType is LINEAR_STRETCHING
    1107       else if( styleType == COLOR_PALETTE )
     1168      else if( styleType == CLASSIFY )
    11081169        {
    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)
     1170        if(iBand==0)
    11271171          {
    1128           double delta = max - min;
    1129           double interval = delta / 10;
    1130           double cstep = min;
    1131           for(i=0; i<10; i++)
     1172          if (classifyType == USER)
    11321173            {
    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);
     1174#ifdef DEBUGMS
     1175            fprintf(stderr,"Processing the raster using the style given in %s\n",msRasterStyleOptionsContent);
     1176#endif
     1177
     1178            // Read the mapfile
     1179            FILE *externalMapfile;
     1180            externalMapfile = fopen(msRasterStyleOptionsContent, "rb");
     1181            if (externalMapfile != NULL)
     1182              {
     1183              long lSize;
     1184              char *buffer;
     1185
     1186              fseek( externalMapfile , 0L , SEEK_END);
     1187              lSize = ftell( externalMapfile );
     1188              rewind( externalMapfile );
     1189
     1190              /* allocate memory for entire content */
     1191              buffer = calloc( 1, lSize+1 );
     1192              if( !buffer )
     1193                {
     1194                fprintf(stderr,"Unable to allocate buffer for file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
     1195                classifyType = defaultClassifyType;
     1196                }
     1197              else
     1198                {
     1199
     1200                /* copy the file into the buffer */
     1201                if( 1!=fread( buffer , lSize, 1 , externalMapfile) )
     1202                  {
     1203                  fclose(externalMapfile);
     1204                  free(buffer);
     1205                  fprintf(stderr,"Unable to read entire file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
     1206                  classifyType = defaultClassifyType;
     1207                  }
     1208                else
     1209                  {
     1210
     1211                  /* do your work here, buffer is a string contains the whole text */
     1212
     1213                  fclose(externalMapfile);
     1214                  msUpdateLayerFromString(myLayer, buffer, 0);
     1215                  free(buffer);
     1216                  } // can read entire file
     1217                } // can allocate buffer
     1218              } // file exist
    11581219            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
     1220              {
     1221              fprintf(stderr,"Unable to read file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
     1222              classifyType = defaultClassifyType;
     1223              } // file doesn't exist
     1224
     1225            } // classify type is USER
     1226
     1227          if (classifyType == AUTO)
     1228            {
     1229            // The number of classes is min(delta, maxNbOfClasses)
     1230            double delta = max - min;
     1231            double step = 1.0;
     1232            double lowBound = 1.0 * min;
     1233            unsigned int numberOfClasses = msRasterStyleClassifyAutoMaximumNumberOfClasses;
     1234            if (delta < msRasterStyleClassifyAutoMaximumNumberOfClasses)
     1235              {
     1236              numberOfClasses = (unsigned int) delta + 1;
     1237              }
     1238            else
     1239              {
     1240              step = delta / (1.0 * msRasterStyleClassifyAutoMaximumNumberOfClasses);
     1241              }
     1242#ifdef DEBUGMS
     1243            fprintf(stderr,"Processing the raster using %d classes with values from %.3f with a step of %.3f\n",numberOfClasses, lowBound, step);
     1244#endif
     1245
     1246            for(i=0; i<numberOfClasses; i++)
     1247              {
     1248              /**
     1249               * Create a new class
     1250               */
     1251              if(msGrowLayerClasses(myLayer) == NULL)
     1252                return -1;
     1253              if(initClass((myLayer->CLASS[myLayer->numclasses])) == -1)
     1254                return -1;
     1255              if(msGrowClassStyles(myLayer->CLASS[myLayer->numclasses]) == NULL)
     1256                return -1;
     1257              if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
     1258                return -1;
     1259
     1260              /**
     1261               * Set class name
     1262               */
     1263              char className[7];
     1264              sprintf(className,"class%d",i);
     1265              myLayer->CLASS[myLayer->numclasses]->name=zStrdup(className);
     1266
     1267              /**
     1268               * Set expression
     1269               */
     1270              char expression[1024];
     1271              if(i+1<numberOfClasses)
     1272                sprintf(expression,"([pixel]>=%.3f AND [pixel]<%.3f)",lowBound,lowBound+step);
     1273              else
     1274                sprintf(expression,"([pixel]>=%.3f AND [pixel]<=%.3f)",lowBound,lowBound+step);
     1275              msLoadExpressionString(&myLayer->CLASS[myLayer->numclasses]->expression,expression);
     1276              lowBound += step;
     1277
     1278              /**
     1279               * Set color
     1280               */
     1281              double g = i / (0.5*numberOfClasses) ; // must be in [-1,1]
     1282              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.red=(int)(255*base(g-0.5));
     1283              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.green=(int)(255*base(g));
     1284              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.blue=(int)(255*base(g+0.5));
     1285              myLayer->CLASS[myLayer->numclasses]->numstyles++;
     1286              myLayer->numclasses++;
     1287
     1288              } // next class
     1289
     1290            } // classify type is AUTO
     1291
     1292          } // styleType is CLASSIFY
     1293        } //iBand is 0
    11791294
    11801295      } // If no error with GDAL functions
    … …  
    11941309    if (hasNoData)
    11951310      {
     1311#ifdef DEBUGMS
     1312      fprintf(stderr,"No data detected (%.3f)\n", noDataValue);
     1313#endif
    11961314      offsiteR = (int) noDataValue;
    11971315      offsiteG = (int) noDataValue;
    … …  
    12011319    myLayer->offsite.green  = offsiteG;
    12021320    myLayer->offsite.blue   = offsiteB;
     1321#ifdef DEBUGMS
     1322    fprintf(stderr,"Setting OFFSITE to (%d,%d,%d)\n",offsiteR, offsiteG, offsiteB);
     1323#endif
    12031324
    12041325    /*
    … …  
    12171338   * Check if there is resample option
    12181339   */
    1219   char msResampleOptionInstruction[1024];
     1340  char msResampleOptionDirective[1024];
    12201341  char * msRasterResampleOptionContent = "BILINEAR";
    1221   map* msRasterResamplingOption=getMap(output->content, msRasterResamplingType);
     1342  map* msRasterResamplingOption=getMap(output->content, msRasterResamplingPropertyName);
    12221343  if (msRasterResamplingOption!=NULL)
    12231344    {
    12241345    msRasterResampleOptionContent = msRasterResamplingOption->value;
    12251346    }
    1226   sprintf(msResampleOptionInstruction, "RESAMPLE=%s",msRasterResampleOptionContent);
    1227   msLayerAddProcessing(myLayer, msResampleOptionInstruction);
     1347  sprintf(msResampleOptionDirective, "RESAMPLE=%s",msRasterResampleOptionContent);
     1348  msLayerAddProcessing(myLayer, msResampleOptionDirective);
    12281349
    12291350  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