source: branches/ms-style/zoo-project/zoo-kernel/service_internal_ms.c @ 910

Last change on this file since 910 was 833, checked in by cresson, 8 years ago

DOC: Add some comments
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 48.7 KB
Line 
1/*
2 * Author : Gérald FENOY
3 *
4 *  Copyright 2010-2011 Fondazione Edmund Mach. All rights reserved.
5 *
6 * Permission is hereby granted, free of charge, to any person obtaining a copy
7 * of this software and associated documentation files (the "Software"), to deal
8 * in the Software without restriction, including without limitation the rights
9 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10 * copies of the Software, and to permit persons to whom the Software is
11 * furnished to do so, subject to the following conditions:
12 *
13 * The above copyright notice and this permission notice shall be included in
14 * all copies or substantial portions of the Software.
15 *
16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22 * THE SOFTWARE.
23 */
24
25/**
26 * Cross platform definition of layerObj->class
27 */
28#ifndef WIN32
29#define CLASS class
30#else
31#define CLASS _class
32#endif
33#include "service_internal_ms.h"
34#include "server_internal.h"
35#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 and sub-styles
118 */
119enum MS_RASTER_STYLES{
120  LINEAR_STRETCHING,
121  CLASSIFY
122};
123enum LINEAR_STRETCHING_TYPE{
124  MINMAX,
125  MEANSTD
126};
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}
138double base( double val ) {
139    if ( val <= 0.25 ) return 0;
140    else if ( val <= 0.75 ) return interpolate( val, 0.0, 0.25, 1.0, 0.75 );
141    else if ( val <= 1.25 ) return 1.0;
142    else if ( val <= 1.75 ) return interpolate( val, 1.0, 1.25, 0.0, 1.75 );
143    else return 0.0;
144}
145
146/**
147 * Get a list of configuration keys having a corresponding mandatory ows_*.
148 * Map composed by a main.cfg maps name as key and the corresponding
149 * MapServer Mafile Metadata name to use
150 * see doc from here :
151 *  - http://mapserver.org/ogc/wms_server.html
152 *  - http://mapserver.org/ogc/wfs_server.html
153 *  - http://mapserver.org/ogc/wcs_server.html
154 *
155 * @return a new map containing a table linking a name of a configuration key
156 * to a corresponding mandatory ows_* keyword (ie. "fees" => "ows_fees").
157 */
158map* getCorrespondance(){
159  map* res=createMap("encoding","ows_encoding");
160  addToMap(res,"abstract","ows_abstract");
161  addToMap(res,"title","ows_title");
162  addToMap(res,"keywords","ows_keywordlist");
163  addToMap(res,"fees","ows_fees");
164  addToMap(res,"accessConstraints","ows_accessconstraints");
165  addToMap(res,"providerName","ows_attribution_title");
166  addToMap(res,"providerSite","ows_service_onlineresource");
167  addToMap(res,"individualName","ows_contactperson");
168  addToMap(res,"positionName","ows_contactposition");
169  addToMap(res,"providerName","ows_contactorganization");
170  addToMap(res,"role","ows_role");
171  addToMap(res,"addressType","ows_addresstype");
172  addToMap(res,"addressCity","ows_city");
173  addToMap(res,"addressDeliveryPoint","ows_address");
174  addToMap(res,"addressPostalCode","ows_postcode");
175  addToMap(res,"addressAdministrativeArea","ows_stateorprovince");
176  addToMap(res,"addressCountry","ows_country");
177  addToMap(res,"phoneVoice","ows_contactvoicetelephone");
178  addToMap(res,"phoneFacsimile","ows_contactfacsimiletelephone");
179  addToMap(res,"addressElectronicMailAddress","ows_contactelectronicmailaddress");
180  // Missing Madatory Informations
181  addToMap(res,"hoursOfService","ows_hoursofservice");
182  addToMap(res,"contactInstructions","ows_contactinstructions");
183  return res;
184}
185
186/**
187 * Add width and height keys to an output maps containing the maximum width
188 * and height for displaying the full data extent.
189 * Restriction to an image having a size of 640x480 (width * height)
190 *
191 * @param output
192 * @param minx the lower left x coordinate
193 * @param miny the lower left y coordinate
194 * @param maxx the upper right x coordinate
195 * @param maxy the upper right y coordinate
196 */
197void setMapSize(maps* output,double minx,double miny,double maxx,double maxy){
198  double maxWidth=640;
199  double maxHeight=480;
200  double deltaX=maxx-minx;
201  double deltaY=maxy-miny;
202  double qWidth;
203  qWidth=maxWidth/deltaX;
204  double qHeight;
205  qHeight=maxHeight/deltaY;
206#ifdef DEBUGMS
207  fprintf(stderr,"deltaX : %.15f \ndeltaY : %.15f\n",deltaX,deltaY);
208  fprintf(stderr,"qWidth : %.15f \nqHeight : %.15f\n",qWidth,qHeight);
209#endif
210
211  double width=deltaX*qWidth;
212  double height=height=deltaY*qWidth;
213  if(deltaX<deltaY){
214    width=deltaX*qHeight;
215    height=deltaY*qHeight;
216  }
217  if(height<0)
218    height=-height;
219  if(width<0)
220    width=-width;
221  char sWidth[1024];
222  char sHeight[1024];
223  sprintf(sWidth,"%.3f",width);
224  sprintf(sHeight,"%.3f",height);
225#ifdef DEBUGMS
226  fprintf(stderr,"sWidth : %.15f \nsHeight : %.15f\n",sWidth,sHeight);
227#endif
228  if(output!=NULL){
229    addToMap(output->content,"width",sWidth);
230    addToMap(output->content,"height",sHeight);
231  }
232}
233
234/**
235 * Add a Reference key to an output containing the WMFS/WFS/WCS request for
236 * accessing service result
237 *
238 * @param m the conf maps containing the main.cfg settings
239 * @param tmpI the specific output maps to add the Reference key
240 */
241void setReferenceUrl(maps* m,maps* tmpI){
242  outputMapfile(m,tmpI);
243  map *msUrl=getMapFromMaps(m,"main","mapserverAddress");
244  if(msUrl==NULL){
245    errorException (m, _("Unable to find any mapserverAddress defined in the main.cfg file"),
246                    "InternalError", NULL);
247    exit(-1);
248  }
249  map *msOgcVersion=getMapFromMaps(m,"main","msOgcVersion");
250  map *dataPath=getMapFromMaps(m,"main","dataPath");
251  map *sid=getMapFromMaps(m,"lenv","usid");
252  map* format=getMap(tmpI->content,"mimeType");
253  map* rformat=getMap(tmpI->content,"requestedMimeType");
254  map* width=getMap(tmpI->content,"width");
255  map* height=getMap(tmpI->content,"height");
256  map* protoMap=getMap(tmpI->content,"msOgc");
257  map* versionMap=getMap(tmpI->content,"msOgcVersion");
258  char options[3][5][25]={
259    {"WMS","1.3.0","GetMap","layers=%s","wms_extent"},
260    {"WFS","1.1.0","GetFeature","typename=%s","wcs_extent"},
261    {"WCS","1.1.0","GetCoverage","coverage=%s","wcs_extent"}
262  };
263  int proto=0;
264  if(rformat==NULL){
265    rformat=getMap(tmpI->content,"mimeType");
266  }
267  if(strncasecmp(rformat->value,"text/xml",8)==0)
268    proto=1;
269  if(strncasecmp(rformat->value,"image/tiff",10)==0)
270    proto=2;
271  if(protoMap!=NULL){
272    if(strncasecmp(protoMap->value,"WMS",3)==0)
273      proto=0;
274    else{
275      if(strncasecmp(protoMap->value,"WFS",3)==0)
276        proto=1;
277      else 
278        proto=2;
279    }
280  }
281  char *protoVersion=options[proto][1];
282  if(proto==1){
283    if(msOgcVersion!=NULL)
284      protoVersion=msOgcVersion->value;
285    if(versionMap!=NULL)
286      protoVersion=versionMap->value;
287  }
288
289  map* extent=getMap(tmpI->content,options[proto][4]);
290  map* crs=getMap(tmpI->content,"crs");
291  int hasCRS=1;
292  if(crs==NULL){
293    crs=getMapFromMaps(m,"main","crs");
294    if(crs==NULL){
295      crs=createMap("crs","epsg:4326");
296      hasCRS=0;
297    }
298  }
299  char layers[128];
300  sprintf(layers,options[proto][3],tmpI->name);
301
302  char* webService_url=(char*)malloc((strlen(msUrl->value)+strlen(format->value)+strlen(tmpI->name)+strlen(width->value)+strlen(height->value)+strlen(extent->value)+256)*sizeof(char));
303
304  if(proto>0){
305    sprintf(webService_url,
306            "%s?map=%s/%s_%s.map&request=%s&service=%s&version=%s&%s&format=%s&bbox=%s&crs=%s",
307            msUrl->value,
308            dataPath->value,
309            tmpI->name,
310            sid->value,
311            options[proto][2],
312            options[proto][0],
313            protoVersion,
314            layers,
315            rformat->value,
316            extent->value,
317            crs->value
318            );
319  }
320  else{
321    sprintf(webService_url,
322            "%s?map=%s/%s_%s.map&request=%s&service=%s&version=%s&%s&width=%s&height=%s&format=%s&bbox=%s&crs=%s",
323            msUrl->value,
324            dataPath->value,
325            tmpI->name,
326            sid->value,
327            options[proto][2],
328            options[proto][0],
329            protoVersion,
330            layers,
331            width->value,
332            height->value,
333            rformat->value,
334            extent->value,
335            crs->value
336            );
337  }
338  if(hasCRS==0){
339    freeMap(&crs);
340    free(crs);
341  }
342  addToMap(tmpI->content,"Reference",webService_url);
343  free(webService_url);
344}
345
346/**
347 * Set projection for a layer in a MAPFILE using Authority Code and Name if
348 * available or fallback to proj4 definition if available or fallback to
349 * default EPSG:4326
350 *
351 * @param output the output maps
352 * @param m the opened mapObj
353 * @param myLayer the layerObj
354 * @param pszProjection a char* containing the SRS definition in WKT format
355 */
356void setSrsInformations(maps* output,mapObj* m,layerObj* myLayer,
357                        char* pszProjection){
358  OGRSpatialReferenceH  hSRS;
359  map* msSrs=NULL;
360  hSRS = OSRNewSpatialReference(NULL);
361  if( pszProjection!=NULL && strlen(pszProjection)>1){
362    if(OSRImportFromWkt( hSRS, &pszProjection ) == CE_None ){
363      char *proj4Str=NULL;
364      if(OSRGetAuthorityName(hSRS,NULL)!=NULL && 
365         OSRGetAuthorityCode(hSRS,NULL)!=NULL){
366        char tmpSrs[20];
367        sprintf(tmpSrs,"%s:%s",
368                OSRGetAuthorityName(hSRS,NULL),OSRGetAuthorityCode(hSRS,NULL));
369        msLoadProjectionStringEPSG(&m->projection,tmpSrs);
370        msLoadProjectionStringEPSG(&myLayer->projection,tmpSrs);
371       
372        char tmpSrss[256];
373        sprintf(tmpSrss,"EPSG:4326 EPSG:900913 EPSG:3857 %s",tmpSrs);
374       
375        msInsertHashTable(&(m->web.metadata), "ows_srs", tmpSrss);
376        msInsertHashTable(&(myLayer->metadata), "ows_srs", tmpSrss);
377       
378#ifdef DEBUGMS
379        fprintf(stderr,"isGeo %b\n\n",OSRIsGeographic(hSRS)==TRUE);
380#endif
381        if(output!=NULL){
382          if(OSRIsGeographic(hSRS)==TRUE)
383            addToMap(output->content,"crs_isGeographic","true");
384          else
385            addToMap(output->content,"crs_isGeographic","false");
386          addToMap(output->content,"crs",tmpSrs);
387        }
388      }
389      else{
390        OSRExportToProj4(hSRS,&proj4Str);
391        if(proj4Str!=NULL && strlen(proj4Str)>0){
392#ifdef DEBUGMS
393          fprintf(stderr,"PROJ (%s)\n",proj4Str);
394#endif
395          msLoadProjectionString(&(m->projection),proj4Str);
396          msLoadProjectionString(&(myLayer->projection),proj4Str);
397          if(output!=NULL){ 
398            if(OSRIsGeographic(hSRS)==TRUE)
399              addToMap(output->content,"crs_isGeographic","true");
400            else
401              addToMap(output->content,"crs_isGeographic","false");
402          }
403          free(proj4Str);
404        }
405        else{
406          msLoadProjectionStringEPSG(&m->projection,"EPSG:4326");
407          msLoadProjectionStringEPSG(&myLayer->projection,"EPSG:4326");
408          if(output!=NULL){
409            addToMap(output->content,"crs_isGeographic","true");
410          }
411        }
412        if(output!=NULL){
413          addToMap(output->content,"crs","EPSG:4326");
414          addToMap(output->content,"real_extent","true");
415        }
416        msInsertHashTable(&(m->web.metadata),"ows_srs", "EPSG:4326 EPSG:900913 EPSG:3857");
417        msInsertHashTable(&(myLayer->metadata),"ows_srs","EPSG:4326 EPSG:900913 EPSG:3857");
418      }
419    }
420  }
421  else{
422    if(output!=NULL){
423      msSrs=getMap(output->content,"msSrs");
424    }
425    if(msSrs!=NULL){
426      msLoadProjectionStringEPSG(&m->projection,msSrs->value);
427      msLoadProjectionStringEPSG(&myLayer->projection,msSrs->value);
428      char tmpSrs[128];
429      sprintf(tmpSrs,"%s EPSG:4326 EPSG:900913 EPSG:3857",msSrs->value);
430      msInsertHashTable(&(m->web.metadata),"ows_srs",tmpSrs);
431      msInsertHashTable(&(myLayer->metadata),"ows_srs",tmpSrs);
432    }else{
433      msLoadProjectionStringEPSG(&m->projection,"EPSG:4326");
434      msLoadProjectionStringEPSG(&myLayer->projection,"EPSG:4326");
435      msInsertHashTable(&(m->web.metadata),"ows_srs","EPSG:4326 EPSG:900913 EPSG:3857");
436      msInsertHashTable(&(myLayer->metadata),"ows_srs","EPSG:4326 EPSG:900913 EPSG:3857");
437    }
438    if(output!=NULL){
439      addToMap(output->content,"crs",msSrs->value);
440      addToMap(output->content,"crs_isGeographic","true");
441    }
442  }
443
444  OSRDestroySpatialReference( hSRS );
445}
446
447/**
448 * Set the MAPFILE extent, the the ows_extent for the layer, add wms_extent and
449 * wfs_extent to the output maps and call setMapSize.
450 *
451 * @param output the specific output
452 * @param m the mapObj
453 * @param myLayer the layerObj
454 * @param minX the lower left x coordinate
455 * @param minY the lower left y coordinate
456 * @param maxX the upper right x coordinate
457 * @param maxY the upper right y coordinate
458 * @see setMapSize
459 */
460void setMsExtent(maps* output,mapObj* m,layerObj* myLayer,
461                 double minX,double minY,double maxX,double maxY){
462  msMapSetExtent(m,minX,minY,maxX,maxY);
463#ifdef DEBUGMS
464  fprintf(stderr,"Extent %.15f %.15f %.15f %.15f\n",minX,minY,maxX,maxY);
465#endif
466  char tmpExtent[1024];
467  sprintf(tmpExtent,"%.15f %.15f %.15f %.15f",minX,minY,maxX,maxY);
468#ifdef DEBUGMS
469  fprintf(stderr,"Extent %s\n",tmpExtent);
470#endif
471  msInsertHashTable(&(myLayer->metadata), "ows_extent", tmpExtent);
472 
473  if(output!=NULL){
474    map* test=getMap(output->content,"real_extent");
475    if(test!=NULL){
476      pointObj min, max;
477      projectionObj tempSrs;
478      min.x = m->extent.minx;
479      min.y = m->extent.miny;
480      max.x = m->extent.maxx;
481      max.y = m->extent.maxy;
482      char tmpSrsStr[1024];
483      msInitProjection(&tempSrs);
484      msLoadProjectionStringEPSG(&tempSrs,"EPSG:4326");
485
486      msProjectPoint(&(m->projection),&tempSrs,&min);
487      msProjectPoint(&m->projection,&tempSrs,&max);
488     
489      sprintf(tmpExtent,"%.3f,%.3f,%.3f,%.3f",min.y,min.x,max.y,max.x);
490      map* isGeo=getMap(output->content,"crs_isGeographic");
491#ifdef DEBUGMS
492      fprintf(stderr,"isGeo = %s\n",isGeo->value);
493#endif
494      if(isGeo!=NULL && strcasecmp("true",isGeo->value)==0)
495        sprintf(tmpExtent,"%f,%f,%f,%f", minY,minX, maxY, maxX);
496      addToMap(output->content,"wms_extent",tmpExtent);
497      sprintf(tmpSrsStr,"%.3f,%.3f,%.3f,%.3f",min.x,min.y,max.x,max.y);
498      addToMap(output->content,"wcs_extent",tmpExtent);
499    }else{
500      sprintf(tmpExtent,"%f,%f,%f,%f",minX, minY, maxX, maxY);
501      map* isGeo=getMap(output->content,"crs_isGeographic");
502      if(isGeo!=NULL){
503#ifdef DEBUGMS
504        fprintf(stderr,"isGeo = %s\n",isGeo->value);
505#endif
506        if(isGeo!=NULL && strcasecmp("true",isGeo->value)==0)
507          sprintf(tmpExtent,"%f,%f,%f,%f", minY,minX, maxY, maxX);
508      }
509      addToMap(output->content,"wms_extent",tmpExtent); 
510      sprintf(tmpExtent,"%.3f,%.3f,%.3f,%.3f",minX,minY,maxX,maxY);
511      addToMap(output->content,"wcs_extent",tmpExtent);
512    }
513  }
514
515  setMapSize(output,minX,minY,maxX,maxY);
516}
517
518/**
519 * Try to open a vector output and define the corresponding layer in the MAPFILE
520 *
521 * @param conf the conf maps containing the main.cfg settings
522 * @param output the specific output maps
523 * @param m the mapObj
524 */
525int tryOgr(maps* conf,maps* output,mapObj* m){
526
527  map* tmpMap=getMap(output->content,"storage");
528  char *pszDataSource=tmpMap->value;
529
530  /**
531   * Try to open the DataSource using OGR
532   */
533  OGRRegisterAll();
534  /**
535   * Try to load the file as ZIP
536   */
537
538  OGRDataSourceH poDS1 = NULL;
539  OGRSFDriverH *poDriver1 = NULL;
540  char *dsName=(char*)malloc((8+strlen(pszDataSource)+1)*sizeof(char));
541  char *odsName=zStrdup(pszDataSource);
542  char *sdsName=zStrdup(pszDataSource);
543  char *demo=".data";
544  sdsName[strlen(sdsName)-(strlen(demo)-1)]='d';
545  sdsName[strlen(sdsName)-(strlen(demo)-2)]='i';
546  sdsName[strlen(sdsName)-(strlen(demo)-3)]='r';
547  sdsName[strlen(sdsName)-(strlen(demo)-4)]=0;
548
549  odsName[strlen(odsName)-(strlen(demo)-1)]='z';
550  odsName[strlen(odsName)-(strlen(demo)-2)]='i';
551  odsName[strlen(odsName)-(strlen(demo)-3)]='p';
552  odsName[strlen(odsName)-(strlen(demo)-4)]=0;
553  sprintf(dsName,"/vsizip/%s",odsName);
554
555#ifdef DEBUGMS
556  fprintf(stderr,"Try loading %s, %s, %s\n",dsName,odsName,dsName);
557#endif
558
559  FILE* file = fopen(pszDataSource, "rb");
560  FILE* fileZ = fopen(odsName, "wb");
561  free(odsName);
562  fseek(file, 0, SEEK_END);
563  unsigned long fileLen=ftell(file);
564  fseek(file, 0, SEEK_SET);
565  char *buffer=(char *)malloc(fileLen+1);
566  fread(buffer, fileLen, 1, file);
567  fwrite(buffer,fileLen, 1, fileZ);
568  fclose(file);
569  fclose(fileZ);
570  free(buffer);
571#ifdef DEBUGMS
572  fprintf(stderr,"Try loading %s",dsName);
573#endif
574  poDS1 = OGROpen( dsName, FALSE, poDriver1 );
575  if( poDS1 == NULL ){
576    fprintf(stderr,"Unable to access the DataSource as ZIP File\n");
577    setMapInMaps(conf,"lenv","message","Unable to open datasource in read only mode");
578    OGR_DS_Destroy(poDS1);
579  }else{
580#ifdef DEBUGMS
581    fprintf(stderr,"The DataSource is a  ZIP File\n");
582#endif
583    char** demo=VSIReadDir(dsName);
584    int i=0;
585    zMkdir(sdsName
586#ifndef WIN32
587                ,S_IRWXU | S_IRGRP | S_IXGRP | S_IROTH | S_IXOTH
588#endif
589                );
590    while(demo[i]!=NULL){
591#ifdef DEBUGMS
592      fprintf(stderr,"ZIP File content : %s\n",demo[i]);
593#endif
594      char *tmpDs=(char*)malloc((strlen(dsName)+strlen(demo[i])+2)*sizeof(char));
595      sprintf(tmpDs,"%s/%s",dsName,demo[i]);
596      fprintf(stderr,"read : %s\n",tmpDs);
597     
598      VSILFILE* vsif=VSIFOpenL(tmpDs,"rb");
599#ifdef DEBUGMS
600      fprintf(stderr,"open : %s\n",tmpDs);
601#endif
602      VSIFSeekL(vsif,0,SEEK_END);
603      vsi_l_offset size=VSIFTellL(vsif);
604#ifdef DEBUGMS
605      fprintf(stderr,"size : %d\n",size);
606#endif
607      VSIFSeekL(vsif,0,SEEK_SET);
608      char *vsifcontent=(char*) malloc(((int)size+1)*sizeof(char));
609      VSIFReadL(vsifcontent,1,(size_t)size,vsif);
610      char *fpath=(char*) malloc((strlen(sdsName)+strlen(demo[1])+2)*sizeof(char));
611      sprintf(fpath,"%s/%s",sdsName,demo[i]);
612      int f=zOpen(fpath,O_WRONLY|O_CREAT,S_IRUSR|S_IWUSR|S_IRGRP|S_IWGRP|S_IROTH|S_IWOTH);
613      zWrite(f,vsifcontent,(int)size);
614      close(f);
615      chmod(fpath,S_IRWXU | S_IRGRP | S_IXGRP | S_IROTH);
616      char* tmpP=strstr(fpath,".shp");
617      if(tmpP==NULL)
618        tmpP=strstr(fpath,".SHP");
619      if(tmpP!=NULL){
620#ifdef DEBUGMS
621        fprintf(stderr,"*** DEBUG %s\n",strstr(tmpP,"."));
622#endif
623        if( strcmp(tmpP,".shp")==0 || strcmp(tmpP,".SHP")==0 ){
624          tmpMap=getMap(output->content,"storage");
625          free(tmpMap->value);
626          tmpMap->value=(char*) malloc((strlen(fpath)+1)*sizeof(char));
627          sprintf(tmpMap->value,"%s",fpath);
628          pszDataSource=tmpMap->value;
629#ifdef DEBUGMS
630          fprintf(stderr,"*** DEBUG %s\n",pszDataSource);
631#endif
632        }
633      }
634      VSIFCloseL(vsif);
635      i++;
636    }
637  }
638  free(sdsName);
639  free(dsName);
640
641  OGRDataSourceH poDS = NULL;
642  OGRSFDriverH *poDriver = NULL;
643  poDS = OGROpen( pszDataSource, FALSE, poDriver );
644  if( poDS == NULL ){
645#ifdef DEBUGMS
646    fprintf(stderr,"Unable to access the DataSource %s\n",pszDataSource);
647#endif
648    setMapInMaps(conf,"lenv","message","Unable to open datasource in read only mode");
649    OGR_DS_Destroy(poDS);
650    OGRCleanupAll();
651#ifdef DEBUGMS
652    fprintf(stderr,"Unable to access the DataSource, exit! \n"); 
653#endif
654    return -1;
655  }
656
657  int iLayer = 0;
658  for( iLayer=0; iLayer < OGR_DS_GetLayerCount(poDS); iLayer++ ){
659    OGRLayerH poLayer = OGR_DS_GetLayer(poDS,iLayer);
660
661    if( poLayer == NULL ){
662#ifdef DEBUGMS
663      fprintf(stderr,"Unable to access the DataSource Layer \n");
664#endif
665      setMapInMaps(conf,"lenv","message","Unable to open datasource in read only mode");
666      return -1;
667    }
668
669    /**
670     * Add a new layer set name, data
671     */
672    if(msGrowMapLayers(m)==NULL){
673      return -1;
674    }
675    if(initLayer((m->layers[m->numlayers]), m) == -1){
676      return -1;
677    }
678
679    layerObj* myLayer=m->layers[m->numlayers];
680#ifdef DEBUGMS
681    dumpMaps(output);
682#endif
683    myLayer->name = zStrdup(output->name);
684    myLayer->tileitem=NULL;
685    myLayer->data = zStrdup(OGR_L_GetName(poLayer));
686    myLayer->connection = zStrdup(pszDataSource);
687    myLayer->index = m->numlayers;
688    myLayer->dump = MS_TRUE;
689    myLayer->status = MS_ON;
690    msConnectLayer(myLayer,MS_OGR,pszDataSource);
691
692    /**
693     * Detect the Geometry Type or use Polygon
694     */
695    if(OGR_L_GetGeomType(poLayer) != wkbUnknown){
696      switch(OGR_L_GetGeomType(poLayer)){
697      case wkbPoint:
698      case wkbMultiPoint:
699      case wkbPoint25D:
700      case wkbMultiPoint25D:
701#ifdef DEBUGMS
702        fprintf(stderr,"POINT DataSource Layer \n");
703#endif
704        myLayer->type = MS_LAYER_POINT;
705        break;
706      case wkbLineString :
707      case wkbMultiLineString :
708      case wkbLineString25D:
709      case wkbMultiLineString25D:
710#ifdef DEBUGMS
711        fprintf(stderr,"LINE DataSource Layer \n");
712#endif
713        myLayer->type = MS_LAYER_LINE;
714        break;
715      case wkbPolygon:
716      case wkbMultiPolygon:
717      case wkbPolygon25D:
718      case wkbMultiPolygon25D:
719#ifdef DEBUGMS
720        fprintf(stderr,"POLYGON DataSource Layer \n");
721#endif
722        myLayer->type = MS_LAYER_POLYGON;
723        break;
724      default:
725        myLayer->type = MS_LAYER_POLYGON;
726        break;
727      }
728    }else
729      myLayer->type = MS_LAYER_POLYGON;
730
731    /**
732     * Detect spatial reference or use WGS84
733     */
734    OGRSpatialReferenceH srs=OGR_L_GetSpatialRef(poLayer);
735    if(srs!=NULL){
736      char *wkt=NULL;
737      OSRExportToWkt(srs,&wkt);
738      setSrsInformations(output,m,myLayer,wkt);
739      free(wkt);
740    }
741    else{
742      addToMap(output->content,"crs","EPSG:4326");
743      addToMap(output->content,"crs_isGeographic","true");
744      msLoadProjectionStringEPSG(&m->projection,"EPSG:4326");
745      msInsertHashTable(&(m->web.metadata), "ows_srs", "EPSG:4326 EPSG:900913 EPSG:3857");
746      msInsertHashTable(&(myLayer->metadata), "ows_srs", "EPSG:4326 EPSG:900913 EPSG:3857");
747    }
748
749    OGREnvelope oExt;
750    if (OGR_L_GetExtent(poLayer,&oExt, TRUE) == OGRERR_NONE){
751      setMsExtent(output,m,myLayer,oExt.MinX, oExt.MinY, oExt.MaxX, oExt.MaxY);
752    }
753 
754    /**
755     * Detect the FID column or use the first attribute field as FID
756     */
757    char *fid=(char*)OGR_L_GetFIDColumn(poLayer);
758    if(strlen(fid)==0){
759      OGRFeatureDefnH def=OGR_L_GetLayerDefn(poLayer);
760      int fIndex=0;
761      for(fIndex=0;fIndex<OGR_FD_GetFieldCount(def);fIndex++){
762        OGRFieldDefnH fdef=OGR_FD_GetFieldDefn(def,fIndex);
763        fid=(char*)OGR_Fld_GetNameRef(fdef);
764        break;
765      }
766    }
767    msInsertHashTable(&(myLayer->metadata), "gml_featureid", fid);
768    msInsertHashTable(&(myLayer->metadata), "gml_include_items", "all");
769    msInsertHashTable(&(myLayer->metadata), "ows_name", output->name);
770    map* tmpMap=getMap(output->content,"title");
771    if(tmpMap!=NULL)
772      msInsertHashTable(&(myLayer->metadata), "ows_title", tmpMap->value);
773    else
774      msInsertHashTable(&(myLayer->metadata), "ows_title", "Default Title");
775   
776    if(msGrowLayerClasses(myLayer) == NULL)
777      return -1;
778    if(initClass((myLayer->CLASS[myLayer->numclasses])) == -1)
779      return -1;
780    if(msGrowClassStyles(myLayer->CLASS[myLayer->numclasses]) == NULL)
781      return -1;
782    if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
783      return -1;
784
785    /**
786     * Apply msStyle else fallback to the default style
787     */
788    tmpMap=getMap(output->content,"msStyle");
789    if(tmpMap!=NULL)
790      msUpdateStyleFromString(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles],tmpMap->value,0);
791    else{
792      /**
793       * Set style
794       */
795      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.red=125;
796      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.green=125;
797      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.blue=255;
798      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->outlinecolor.red=80;
799      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->outlinecolor.green=80;
800      myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->outlinecolor.blue=80;
801     
802      /**
803       * Set specific style depending on type
804       */
805      if(myLayer->type == MS_LAYER_POLYGON)
806        myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->width=3;
807      if(myLayer->type == MS_LAYER_LINE){
808        myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->width=3;
809        myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->outlinewidth=1.5;
810      }
811      if(myLayer->type == MS_LAYER_POINT){
812        myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->symbol=1;
813        myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->size=15;
814      }
815     
816    }
817    myLayer->CLASS[myLayer->numclasses]->numstyles++;
818    myLayer->numclasses++;
819   
820    m->layerorder[m->numlayers] = m->numlayers;
821    m->numlayers++;
822
823  }
824
825  OGR_DS_Destroy(poDS);
826  OGRCleanupAll();
827
828  return 1;
829}
830
831/**
832 * Try to open a raster output and define the corresponding layer in the MAPFILE
833 *
834 * @param conf the conf maps containing the main.cfg settings
835 * @param output the specific output maps
836 * @param m the mapObj
837 */
838int tryGdal(maps* conf,maps* output,mapObj* m){
839
840  /*
841   * Detect the raster style
842   */
843
844  // Properties names
845  const char * msRasterResamplingPropertyName       = "msRasterResample";       /* can be NEAREST/AVERAGE/BILINEAR */
846  const char * msRasterStylePropertyName            = "msRasterStyle";          /* can be linearStretching/classify */
847  const char * msRasterStyleOptionsPropertyName     = "msRasterStyleOptions";   /* can be the path to a mapfile */
848
849  // Allowed properties values
850  const char * msRasterStyleLinearStretchingPropertyValue       = "linearStretching";   /* RasterStyle */
851  const char * msRasterStyleColorPalettePropertyValue           = "classify";
852  const char * msRasterStyleLinearStretchingMinMaxPropertyName  = "minMax";             /* Linear stretching */
853  const char * msRasterStyleLinearStretchingMeanStdPropertyName = "meanStd";
854
855  // Maximum number of classes (classify mode)
856  const unsigned int msRasterStyleClassifyAutoMaximumNumberOfClasses = 256;
857
858  // Default raster style
859  int defaultStyleType = LINEAR_STRETCHING;
860  int defaultLinearstretchingType = MEANSTD;
861  int defaultClassifyType = AUTO;
862
863  // Check if there is a defined raster style type
864  int styleType = defaultStyleType;
865  int linearStretchingType = defaultLinearstretchingType;
866  int classifyType = defaultClassifyType;
867  map* msRasterStyle=getMap(output->content, msRasterStylePropertyName);
868  char * msRasterStyleOptionsContent = "";
869  char * msRasterStyleFileContent = "";
870  if(msRasterStyle!=NULL)
871    {
872#ifdef DEBUGMS
873    fprintf(stderr,"msRasterStyle=%s\n", msRasterStyle->value);
874#endif
875    // Check if there is options attached
876    map* msRasterStyleOptions=getMap(output->content, msRasterStyleOptionsPropertyName);
877    if (msRasterStyleOptions!=NULL)
878      {
879      msRasterStyleOptionsContent = msRasterStyleOptions->value;
880#ifdef DEBUGMS
881      fprintf(stderr,"msRasterStyleOptions=%s\n", msRasterStyleOptionsContent);
882#endif
883
884      }
885
886    if (strncasecmp(msRasterStyle->value, msRasterStyleLinearStretchingPropertyValue,
887        strlen(msRasterStyleLinearStretchingPropertyValue))==0)
888      {
889#ifdef DEBUGMS
890      fprintf(stderr,"The raster style is linear stretching\n");
891#endif
892      styleType = LINEAR_STRETCHING;
893
894      if (strlen(msRasterStyleOptionsContent)>0)
895        {
896        if (strncasecmp(msRasterStyleOptionsContent, msRasterStyleLinearStretchingMinMaxPropertyName,
897            strlen(msRasterStyleLinearStretchingMinMaxPropertyName))==0)
898          {
899          linearStretchingType = MINMAX;
900#ifdef DEBUGMS
901          fprintf(stderr,"The raster style linear stretching option is minmax\n");
902#endif
903          }
904        else if (strncasecmp(msRasterStyleOptionsContent, msRasterStyleLinearStretchingMeanStdPropertyName,
905            strlen(msRasterStyleLinearStretchingMeanStdPropertyName))==0)
906          {
907          linearStretchingType = MEANSTD;
908#ifdef DEBUGMS
909          fprintf(stderr,"The raster style linear stretching option is meanstd\n");
910#endif
911          }
912        else
913          {
914          fprintf(stderr,"Unknown raster style linear stretching method: %s\n", msRasterStyleOptionsContent);
915          }
916        } // raster style options (for linear stretching) are not empty
917      else
918        {
919        fprintf(stderr,"Raster style options for linear stretching are empty. Using default.\n");
920        }
921      } // raster style is linear stretching
922
923    else if (strncasecmp(msRasterStyle->value, msRasterStyleColorPalettePropertyValue,
924        strlen(msRasterStyleColorPalettePropertyValue))==0)
925      {
926#ifdef DEBUGMS
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
946    else
947      {
948      fprintf(stderr,"Unknown raster style: %s. Using default.", msRasterStyle->value);
949      }
950
951    } // raster style is not null
952
953  /*
954   * This is just for the backward compatibility and will be deprecated
955   * TODO: mark this as deprecated
956   */
957  {
958  map* test=getMap(output->content,"msClassify");
959  if(test!=NULL && strncasecmp(test->value,"true",4)==0)
960    {
961    styleType = CLASSIFY;
962    classifyType = AUTO;
963    }
964  }
965
966#ifdef DEBUGMS
967  fprintf(stderr,"Styling options:\n\tRasterStyle=%i\n\tLinearStretching options=%i\n\tClassify options=%i\n",
968      styleType,linearStretchingType, classifyType);
969#endif
970
971  /*
972   * Get storage
973   */
974  map* tmpMap=getMap(output->content,"storage");
975  char *pszFilename=tmpMap->value;
976  GDALDatasetH hDataset;
977  GDALRasterBandH hBand;
978  double adfGeoTransform[6];
979  int i, iBand;
980
981  /**
982   * Try to open the DataSource using GDAL
983   */
984  GDALAllRegister();
985  hDataset = GDALOpen( pszFilename, GA_ReadOnly );
986  if( hDataset == NULL ){
987#ifdef DEBUGMS
988    fprintf(stderr,"Unable to access the DataSource %s \n",pszFilename);
989#endif
990    setMapInMaps(conf,"lenv","message","gdalinfo failed - unable to open");
991    GDALDestroyDriverManager();
992    return -1;
993  }
994#ifdef DEBUGMS
995  fprintf(stderr,"Accessing the DataSource %s %d\n",pszFilename,__LINE__);
996#endif
997
998  /**
999   * Add a new layer set name, data
1000   */
1001  if(msGrowMapLayers(m)==NULL){
1002    return -1;
1003  }
1004  if(initLayer((m->layers[m->numlayers]), m) == -1){
1005    return -1;
1006  }
1007  m->layers[m->numlayers]->index=m->numlayers;
1008
1009  layerObj* myLayer=m->layers[m->numlayers];
1010  myLayer->name = zStrdup(output->name);
1011  myLayer->tileitem=NULL;
1012  myLayer->data = zStrdup(pszFilename);
1013  myLayer->index = m->numlayers;
1014  myLayer->dump = MS_TRUE;
1015  myLayer->status = MS_ON;
1016  myLayer->type = MS_LAYER_RASTER;
1017
1018  char *title=output->name;
1019  tmpMap=getMap(output->content,"title");
1020  if(tmpMap!=NULL)
1021    title=tmpMap->value;
1022  char *abstract=output->name;
1023  tmpMap=getMap(output->content,"abstract");
1024  if(tmpMap!=NULL)
1025    abstract=tmpMap->value;
1026
1027  msInsertHashTable(&(myLayer->metadata), "ows_label", title);
1028  msInsertHashTable(&(myLayer->metadata), "ows_title", title);
1029  msInsertHashTable(&(myLayer->metadata), "ows_abstract", abstract);
1030  msInsertHashTable(&(myLayer->metadata), "ows_rangeset_name", output->name);
1031  msInsertHashTable(&(myLayer->metadata), "ows_rangeset_label", title);
1032
1033  /**
1034   * Set Map Size to the raster size
1035   */
1036  m->width=GDALGetRasterXSize( hDataset );
1037  m->height=GDALGetRasterYSize( hDataset );
1038
1039  /**
1040   * Set projection using Authority Code and Name if available or fallback to
1041   * proj4 definition if available or fallback to default EPSG:4326
1042   */
1043  const char *tRef=GDALGetProjectionRef( hDataset );
1044  if( tRef != NULL && strlen(tRef)>0 ){
1045    char *pszProjection;
1046    pszProjection = (char *) GDALGetProjectionRef( hDataset );
1047#ifdef DEBUGMS
1048    fprintf(stderr,"Accessing the DataSource %s\n",GDALGetProjectionRef( hDataset ));
1049#endif
1050    setSrsInformations(output,m,myLayer,pszProjection);
1051  }else{
1052    fprintf(stderr,"NO SRS FOUND ! %s\n",GDALGetProjectionRef( hDataset ));   
1053  }
1054
1055  /**
1056   * Set extent
1057   */
1058  if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ){
1059    if( adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0 ){
1060
1061      double minX = adfGeoTransform[0]
1062                                    + adfGeoTransform[2] * GDALGetRasterYSize(hDataset);
1063      double minY = adfGeoTransform[3]
1064                                    + adfGeoTransform[5] * GDALGetRasterYSize(hDataset);
1065
1066      double maxX = adfGeoTransform[0]
1067                                    + adfGeoTransform[1] * GDALGetRasterXSize(hDataset);
1068      double maxY = adfGeoTransform[3]
1069                                    + adfGeoTransform[4] * GDALGetRasterXSize(hDataset);
1070
1071      setMsExtent(output,m,myLayer,minX,minY,maxX,maxY);
1072
1073    }
1074  }
1075
1076  /**
1077   * Extract information about available bands to set the bandcount and the
1078   * processing directive
1079   */
1080  char nBands[3];
1081  int nBandsI=GDALGetRasterCount( hDataset );
1082  sprintf(nBands,"%d",GDALGetRasterCount( hDataset ));
1083  msInsertHashTable(&(myLayer->metadata), "ows_bandcount", nBands);
1084  if(nBandsI>=3)
1085    msLayerAddProcessing(myLayer,"BANDS=1,2,3");
1086  else if(nBandsI>=2)
1087    msLayerAddProcessing(myLayer,"BANDS=1,2");
1088  else
1089    msLayerAddProcessing(myLayer,"BANDS=1");
1090
1091  /**
1092   * Name available Bands
1093   */
1094  char lBands[7];
1095  char *nameBands=NULL;
1096  for( iBand = 0; iBand < nBandsI; iBand++ ){
1097    sprintf(lBands,"Band%d",iBand+1);
1098    if(nameBands==NULL){
1099      nameBands=(char*)malloc((strlen(lBands)+1)*sizeof(char));
1100      sprintf(nameBands,"%s",lBands);
1101    }else{
1102      if(iBand<4){
1103        char *tmpS=zStrdup(nameBands);
1104        nameBands=(char*)realloc(nameBands,(strlen(nameBands)+strlen(lBands)+1)*sizeof(char));
1105        sprintf(nameBands,"%s %s",tmpS,lBands);
1106        free(tmpS);
1107      }
1108    }
1109  }
1110  msInsertHashTable(&(myLayer->metadata), "ows_bandnames", nameBands);
1111
1112  /**
1113   * Loops over metadata information to setup specific information
1114   */
1115  for( iBand = 0; iBand < nBandsI; iBand++ )
1116    {
1117
1118    // Compute statistics of the current band
1119    hBand = GDALGetRasterBand( hDataset, iBand+1 );
1120    CPLErrorReset();
1121    double min = 0.0;
1122    double max = 0.0;
1123    double mean = 0.0;
1124    double std = 0.0;
1125    GDALComputeRasterStatistics  (hBand, 1, &min, &max, &mean, &std, NULL, NULL);
1126#ifdef DEBUGMS
1127      fprintf(stderr,"Computed raster stats for band %i: min=%.3f max=%.3f mean=%.3f std=%.3f\n",
1128          iBand, min, max, mean, std);
1129#endif
1130    char bandIdentifier[21];
1131    sprintf(bandIdentifier,"Band%d",iBand+1);
1132    if (CPLGetLastErrorType() == CE_None)
1133      {
1134      char bandInterval[100];
1135      sprintf(bandInterval,"%.3f %.3f",min,max);
1136      char bandIntervalIdentifier[21];
1137      sprintf(bandIntervalIdentifier,"%s_interval",bandIdentifier);
1138      msInsertHashTable(&(myLayer->metadata), bandIntervalIdentifier, bandInterval);
1139
1140      // Apply the raster style
1141      if(styleType == LINEAR_STRETCHING)
1142        {
1143
1144        char msProcessingDirective[1024];
1145        double low = 0.0;
1146        double hi = 1.0;
1147
1148        char s1[MAX_NUMBER_STRING_SIZE];
1149        char s2[MAX_NUMBER_STRING_SIZE];
1150
1151        int bn = iBand+1;
1152        if(linearStretchingType==MINMAX)
1153          {
1154          low = min;
1155          hi = max;
1156          }
1157        else if (linearStretchingType==MEANSTD)
1158          {
1159          low = mean - 2*std;
1160          hi = mean + 2*std;
1161          }
1162#ifdef DEBUGMS
1163      fprintf(stderr,"Processing the raster using a stretch btw %.3f and %.3f\n", low, hi);
1164#endif
1165        sprintf(msProcessingDirective, "SCALE_%d=%s,%s", bn, dtoa(s1,low), dtoa(s2,hi));
1166        msLayerAddProcessing(myLayer,msProcessingDirective);
1167
1168        } // styleType is LINEAR_STRETCHING
1169
1170      else if( styleType == CLASSIFY )
1171        {
1172        if(iBand==0) // classify only the first band
1173          {
1174          if (classifyType == USER)
1175            {
1176#ifdef DEBUGMS
1177            fprintf(stderr,"Processing the raster using the style given in %s\n",msRasterStyleOptionsContent);
1178#endif
1179
1180            // Read the mapfile
1181            FILE *externalMapfile;
1182            externalMapfile = fopen(msRasterStyleOptionsContent, "rb");
1183            if (externalMapfile != NULL)
1184              {
1185              long lSize;
1186              char *buffer;
1187
1188              fseek( externalMapfile , 0L , SEEK_END);
1189              lSize = ftell( externalMapfile );
1190              rewind( externalMapfile );
1191
1192              /* allocate memory for entire content */
1193              buffer = calloc( 1, lSize+1 );
1194              if( !buffer )
1195                {
1196                fprintf(stderr,"Unable to allocate buffer for file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
1197                classifyType = defaultClassifyType;
1198                }
1199              else
1200                {
1201
1202                /* copy the file into the buffer */
1203                if( 1!=fread( buffer , lSize, 1 , externalMapfile) )
1204                  {
1205                  fclose(externalMapfile);
1206                  free(buffer);
1207                  fprintf(stderr,"Unable to read entire file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
1208                  classifyType = defaultClassifyType;
1209                  }
1210                else
1211                  {
1212
1213                  /* do your work here, buffer is a string contains the whole text */
1214
1215                  fclose(externalMapfile);
1216                  msUpdateLayerFromString(myLayer, buffer, 0);
1217                  free(buffer);
1218                  } // can read entire file
1219                } // can allocate buffer
1220              } // file exist
1221            else
1222              {
1223              fprintf(stderr,"Unable to read file %s. Switching to default classes style.\n",msRasterStyleOptionsContent);
1224              classifyType = defaultClassifyType;
1225              } // file doesn't exist
1226
1227            } // classify type is USER
1228
1229          if (classifyType == AUTO)
1230            {
1231            // The number of classes is min(delta, maxNbOfClasses)
1232            double delta = max - min;
1233            double step = 1.0;
1234            double lowBound = 1.0 * min;
1235            unsigned int numberOfClasses = msRasterStyleClassifyAutoMaximumNumberOfClasses;
1236            if (delta < msRasterStyleClassifyAutoMaximumNumberOfClasses)
1237              {
1238              numberOfClasses = (unsigned int) delta + 1;
1239              }
1240            else
1241              {
1242              step = delta / (1.0 * msRasterStyleClassifyAutoMaximumNumberOfClasses);
1243              }
1244#ifdef DEBUGMS
1245            fprintf(stderr,"Processing the raster using %d classes with values from %.3f with a step of %.3f\n",numberOfClasses, lowBound, step);
1246#endif
1247
1248            for(i=0; i<numberOfClasses; i++)
1249              {
1250              /**
1251               * Create a new class
1252               */
1253              if(msGrowLayerClasses(myLayer) == NULL)
1254                return -1;
1255              if(initClass((myLayer->CLASS[myLayer->numclasses])) == -1)
1256                return -1;
1257              if(msGrowClassStyles(myLayer->CLASS[myLayer->numclasses]) == NULL)
1258                return -1;
1259              if(initStyle(myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]) == -1)
1260                return -1;
1261
1262              /**
1263               * Set class name
1264               */
1265              char className[7];
1266              sprintf(className,"class%d",i);
1267              myLayer->CLASS[myLayer->numclasses]->name=zStrdup(className);
1268
1269              /**
1270               * Set expression
1271               */
1272              char expression[1024];
1273              if(i+1<numberOfClasses)
1274                sprintf(expression,"([pixel]>=%.3f AND [pixel]<%.3f)",lowBound,lowBound+step);
1275              else
1276                sprintf(expression,"([pixel]>=%.3f AND [pixel]<=%.3f)",lowBound,lowBound+step);
1277              msLoadExpressionString(&myLayer->CLASS[myLayer->numclasses]->expression,expression);
1278              lowBound += step;
1279
1280              /**
1281               * Set color
1282               */
1283              double g = i / (0.5*numberOfClasses) ; // must be in [-1,1]
1284              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.red=(int)(255*base(g-0.5));
1285              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.green=(int)(255*base(g));
1286              myLayer->CLASS[myLayer->numclasses]->styles[myLayer->CLASS[myLayer->numclasses]->numstyles]->color.blue=(int)(255*base(g+0.5));
1287              myLayer->CLASS[myLayer->numclasses]->numstyles++;
1288              myLayer->numclasses++;
1289
1290              } // next class
1291
1292            } // classify type is AUTO
1293
1294          } // styleType is CLASSIFY
1295
1296        } //iBand is 0
1297
1298      } // If no error with GDAL functions
1299    else
1300      {
1301      fprintf(stderr,"Unable to compute raster statistics!\n");
1302      }
1303
1304    /*
1305     * Set offsite
1306     */
1307    int offsiteR = 0;
1308    int offsiteG = 0;
1309    int offsiteB = 0;
1310    int hasNoData = 0;
1311    double noDataValue = GDALGetRasterNoDataValue(hBand, &hasNoData);
1312    if (hasNoData)
1313      {
1314#ifdef DEBUGMS
1315      fprintf(stderr,"No data detected (%.3f)\n", noDataValue);
1316#endif
1317      offsiteR = (int) noDataValue;
1318      offsiteG = (int) noDataValue;
1319      offsiteB = (int) noDataValue;
1320      }
1321    myLayer->offsite.red    = offsiteR;
1322    myLayer->offsite.green  = offsiteG;
1323    myLayer->offsite.blue   = offsiteB;
1324#ifdef DEBUGMS
1325    fprintf(stderr,"Setting OFFSITE to (%d,%d,%d)\n",offsiteR, offsiteG, offsiteB);
1326#endif
1327
1328    /*
1329     * Insert units
1330     */
1331    if( strlen(GDALGetRasterUnitType(hBand)) > 0 )
1332      {
1333      char tmpU[21];
1334      sprintf(tmpU,"%s_band_uom",bandIdentifier);
1335      msInsertHashTable(&(myLayer->metadata), tmpU, GDALGetRasterUnitType(hBand));
1336      }
1337
1338    } // next band
1339
1340  /*
1341   * Check if there is resample option
1342   */
1343  char msResampleOptionDirective[1024];
1344  char * msRasterResampleOptionContent = "BILINEAR";
1345  map* msRasterResamplingOption=getMap(output->content, msRasterResamplingPropertyName);
1346  if (msRasterResamplingOption!=NULL)
1347    {
1348    msRasterResampleOptionContent = msRasterResamplingOption->value;
1349    }
1350  sprintf(msResampleOptionDirective, "RESAMPLE=%s",msRasterResampleOptionContent);
1351  msLayerAddProcessing(myLayer, msResampleOptionDirective);
1352
1353  m->layerorder[m->numlayers] = m->numlayers;
1354  m->numlayers++;
1355  GDALClose( hDataset );
1356  GDALDestroyDriverManager();
1357  CPLCleanupTLS();
1358  return 1;
1359}
1360
1361/**
1362 * Create a MapFile for WMS, WFS or WCS Service output
1363 *
1364 * @param conf the conf maps containing the main.cfg settings
1365 * @param outputs a specific output maps
1366 */
1367void outputMapfile(maps* conf,maps* outputs){
1368
1369  /**
1370   * First store the value on disk
1371   */
1372  map* mime=getMap(outputs->content,"mimeType");
1373  char *ext="data";
1374  if(mime!=NULL)
1375    if(strncasecmp(mime->value,"application/json",16)==0)
1376      ext="json";
1377
1378  map* tmpMap=getMapFromMaps(conf,"main","dataPath");
1379  map* sidMap=getMapFromMaps(conf,"lenv","usid");
1380  char *pszDataSource=(char*)malloc((strlen(tmpMap->value)+strlen(sidMap->value)+strlen(outputs->name)+17)*sizeof(char));
1381  sprintf(pszDataSource,"%s/ZOO_DATA_%s_%s.%s",tmpMap->value,outputs->name,sidMap->value,ext); 
1382  int f=zOpen(pszDataSource,O_WRONLY|O_CREAT,S_IRUSR|S_IWUSR|S_IRGRP|S_IWGRP|S_IROTH|S_IWOTH);
1383  map *gfile=getMap(outputs->content,"generated_file");
1384  if(gfile!=NULL){
1385    readGeneratedFile(conf,outputs->content,gfile->value);         
1386  }
1387  map* sizeMap=getMap(outputs->content,"size");
1388  map* vData=getMap(outputs->content,"value");
1389  if(sizeMap!=NULL){
1390    zWrite(f,vData->value,atoi(sizeMap->value)*sizeof(char));
1391  }
1392  else{
1393    zWrite(f,vData->value,(strlen(vData->value)+1)*sizeof(char));
1394  }
1395  close(f);
1396  addToMap(outputs->content,"storage",pszDataSource);
1397  free(pszDataSource);
1398
1399  /*
1400   * Create an empty map, set name, default size and extent
1401   */
1402  mapObj *myMap=msNewMapObj();
1403  free(myMap->name);
1404  myMap->name=zStrdup("ZOO-Project_WXS_Server");
1405  msMapSetSize(myMap,2048,2048);
1406  msMapSetExtent(myMap,-1,-1,1,1);
1407 
1408  /*
1409   * Set imagepath and imageurl using tmpPath and tmpUrl from main.cfg
1410   */
1411  map *tmp1=getMapFromMaps(conf,"main","tmpPath");
1412  myMap->web.imagepath=zStrdup(tmp1->value);
1413  tmp1=getMapFromMaps(conf,"main","tmpUrl");
1414  myMap->web.imageurl=zStrdup(tmp1->value);
1415 
1416  /*
1417   * Define supported output formats
1418   */
1419  outputFormatObj *o1=msCreateDefaultOutputFormat(NULL,"AGG/PNG","png");
1420  o1->imagemode=MS_IMAGEMODE_RGBA;
1421  o1->transparent=MS_TRUE;
1422  o1->inmapfile=MS_TRUE;
1423  msAppendOutputFormat(myMap,msCloneOutputFormat(o1));
1424  msFreeOutputFormat(o1);
1425
1426#ifdef USE_KML
1427  outputFormatObj *o2=msCreateDefaultOutputFormat(NULL,"KML","kml");
1428  o2->inmapfile=MS_TRUE; 
1429  msAppendOutputFormat(myMap,msCloneOutputFormat(o2));
1430  msFreeOutputFormat(o2);
1431#endif
1432
1433  outputFormatObj *o3=msCreateDefaultOutputFormat(NULL,"GDAL/GTiff","tiff");
1434  if(!o3)
1435    fprintf(stderr,"Unable to initialize GDAL driver !\n");
1436  else{
1437    o3->imagemode=MS_IMAGEMODE_BYTE;
1438    o3->inmapfile=MS_TRUE; 
1439    msAppendOutputFormat(myMap,msCloneOutputFormat(o3));
1440    msFreeOutputFormat(o3);
1441  }
1442
1443  outputFormatObj *o4=msCreateDefaultOutputFormat(NULL,"GDAL/AAIGRID","grd");
1444  if(!o4)
1445    fprintf(stderr,"Unable to initialize GDAL driver !\n");
1446  else{
1447    o4->imagemode=MS_IMAGEMODE_INT16;
1448    o4->inmapfile=MS_TRUE; 
1449    msAppendOutputFormat(myMap,msCloneOutputFormat(o4));
1450    msFreeOutputFormat(o4);
1451  }
1452
1453#ifdef USE_CAIRO
1454  outputFormatObj *o5=msCreateDefaultOutputFormat(NULL,"CAIRO/PNG","cairopng");
1455  if(!o5)
1456    fprintf(stderr,"Unable to initialize CAIRO driver !\n");
1457  else{
1458    o5->imagemode=MS_IMAGEMODE_RGBA;
1459    o5->transparent=MS_TRUE;
1460    o5->inmapfile=MS_TRUE;
1461    msAppendOutputFormat(myMap,msCloneOutputFormat(o5));
1462    msFreeOutputFormat(o5);
1463  }
1464#endif
1465
1466  /*
1467   * Set default projection to EPSG:4326
1468   */
1469  msLoadProjectionStringEPSG(&myMap->projection,"EPSG:4326");
1470  myMap->transparent=1;
1471
1472  /**
1473   * Set metadata extracted from main.cfg file maps
1474   */
1475  maps* cursor=conf;
1476  map* correspondance=getCorrespondance();
1477  while(cursor!=NULL){
1478    map* _cursor=cursor->content;
1479    map* vMap;
1480    while(_cursor!=NULL){
1481      if((vMap=getMap(correspondance,_cursor->name))!=NULL){
1482        if (msInsertHashTable(&(myMap->web.metadata), vMap->value, _cursor->value) == NULL){
1483#ifdef DEBUGMS
1484          fprintf(stderr,"Unable to add metadata");
1485#endif
1486          return;
1487        }
1488      }
1489      _cursor=_cursor->next;
1490    }
1491    cursor=cursor->next;
1492  }
1493  freeMap(&correspondance);
1494  free(correspondance);
1495
1496  /*
1497   * Set mapserver PROJ_LIB/GDAL_DATA or any other config parameter from
1498   * the main.cfg [mapserver] section if any
1499   */
1500  maps *tmp3=getMaps(conf,"mapserver");
1501  if(tmp3!=NULL){
1502    map* tmp4=tmp3->content;
1503    while(tmp4!=NULL){
1504      msSetConfigOption(myMap,tmp4->name,tmp4->value);
1505      tmp4=tmp4->next;
1506    }
1507  }
1508
1509  /**
1510   * Set a ows_rootlayer_title, 
1511   */
1512  if (msInsertHashTable(&(myMap->web.metadata), "ows_rootlayer_name", "ZOO_Project_Layer") == NULL){
1513#ifdef DEBUGMS
1514    fprintf(stderr,"Unable to add metadata");
1515#endif
1516    return;
1517  }
1518  if (msInsertHashTable(&(myMap->web.metadata), "ows_rootlayer_title", "ZOO_Project_Layer") == NULL){
1519#ifdef DEBUGMS
1520    fprintf(stderr,"Unable to add metadata");
1521#endif
1522    return;
1523  }
1524
1525  /**
1526   * Enable all the WXS requests using ows_enable_request
1527   * see http://mapserver.org/trunk/development/rfc/ms-rfc-67.html
1528   */
1529  if (msInsertHashTable(&(myMap->web.metadata), "ows_enable_request", "*") == NULL){
1530#ifdef DEBUGMS
1531    fprintf(stderr,"Unable to add metadata");
1532#endif
1533    return;
1534  }
1535  msInsertHashTable(&(myMap->web.metadata), "ows_srs", "EPSG:4326");
1536
1537  if(tryOgr(conf,outputs,myMap)<0)
1538    if(tryGdal(conf,outputs,myMap)<0)
1539      return ;
1540
1541  tmp1=getMapFromMaps(conf,"main","dataPath");
1542  char *tmpPath=(char*)malloc((13+strlen(tmp1->value))*sizeof(char));
1543  sprintf(tmpPath,"%s/symbols.sym",tmp1->value);
1544  msInitSymbolSet(&myMap->symbolset);
1545  myMap->symbolset.filename=zStrdup(tmpPath);
1546  free(tmpPath);
1547
1548  map* sid=getMapFromMaps(conf,"lenv","usid");
1549  char *mapPath=
1550    (char*)malloc((7+strlen(sid->value)+strlen(outputs->name)+strlen(tmp1->value))*sizeof(char));
1551  sprintf(mapPath,"%s/%s_%s.map",tmp1->value,outputs->name,sid->value);
1552  msSaveMap(myMap,mapPath);
1553  free(mapPath);
1554  msGDALCleanup();
1555  msFreeMap(myMap);
1556}
1557
Note: See TracBrowser for help on using the repository browser.

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