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

Last change on this file since 832 was 832, checked in by cresson, 7 years ago

ENH: Add custom mapfile styling option

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

Search

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