source: trunk/zoo-project/zoo-services/gdal/profile/service.c @ 975

/* ****************************************************************************
 * $Id$
 *
 * Project:  GdalExtractProfile
 * Purpose:  Extract Profile from a Raster file for an Input Geometry (LINE)
 * Author:   Gérald Fenoy, gerald.fenoy@geolabs.fr
 *
 * ****************************************************************************
 * Copyright (c) 2010-2011, GeoLabs SARL
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#ifdef ZOO_SERVICE
#include "service.h"
#include "service_internal.h"
#endif
#include "gdal.h"
#include "cpl_conv.h"
#include "ogr_api.h"
#if GDAL_VERSION_MAJOR >= 2
#include <gdal_priv.h>
#endif

#ifdef ZOO_SERVICE
extern "C" {
#endif

#ifdef WIN32
__declspec(dllexport)
#endif
#ifdef ZOO_SERVICE
int GdalExtractProfile(maps*& conf,maps*& inputs,maps*& outputs)
#else
int main(int argc,char** argv)
#endif
{
  char *pszFilename;
#ifdef ZOO_SERVICE
  map* tmp=NULL;
  map* tmp1=NULL;
  tmp=getMapFromMaps(conf,"main","dataPath");
  tmp1=getMapFromMaps(inputs,"RasterFile","value");
  pszFilename=(char *)malloc((2+strlen(tmp->value)+strlen(tmp1->value))*sizeof(char));
  sprintf(pszFilename,"%s/%s",tmp->value,tmp1->value);
#else
  pszFilename=argv[1];
#endif
  GDALDatasetH  hDataset;  
  GDALAllRegister();
  OGRRegisterAll();
 
  hDataset = GDALOpen( pszFilename, GA_ReadOnly );
  free(pszFilename);
  if( hDataset != NULL )
    {
      GDALDriverH   hDriver;
      double        adfGeoTransform[6];

      if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
        {


        GDALRasterBandH hBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];
        
        hBand = GDALGetRasterBand( hDataset, 1 );

        adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
        adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
        if( ! (bGotMin && bGotMax) )
          GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );
#ifdef ZOO_SERVICE
        OGRGeometryH geometry=NULL;
        // Verify if there is a cache file for this input, use it if available
        tmp1=getMapFromMaps(inputs,"Geometry","cache_file");
        if(tmp1!=NULL){
#if GDAL_VERSION_MAJOR >= 2
          char* pszDataSource=strdup(tmp1->value);
          GDALDatasetH poDS
            = GDALOpenEx( pszDataSource,
                          GDAL_OF_READONLY | GDAL_OF_VECTOR,
                          NULL, NULL, NULL );
#endif
          for( int iLayer = 0; iLayer < OGR_DS_GetLayerCount(poDS) ;
               iLayer++ )
            {
              OGRLayerH        poLayer = OGR_DS_GetLayer(poDS,iLayer);

              if( poLayer == NULL )
                {
                  fprintf( stderr, "FAILURE: Couldn't fetch advertised layer %d!\n",
                           iLayer );
                  char tmp[1024];
                  sprintf(tmp,"Couldn't fetch advertised layer %d!",iLayer);
                  setMapInMaps(conf,"lenv","message",tmp);
                  return SERVICE_FAILED;
                }
              OGRFeatureH  poFeature;
              int         nFeaturesInTransaction = 0;
              int fCount=0;
              OGR_L_ResetReading(poLayer);
              while( TRUE )
                {
                  poFeature = OGR_L_GetNextFeature(poLayer);
                  if( poFeature == NULL ){
                    break;
                  }
                  geometry=OGR_G_Clone(OGR_F_GetGeometryRef(poFeature));
                  OGR_F_Destroy( poFeature );
                }
              
            }
        }
        else{
          tmp1=getMapFromMaps(inputs,"Geometry","value");
          geometry=OGR_G_CreateGeometryFromJson(tmp1->value);
        }
#else
          OGRGeometryH geometry=OGR_G_CreateGeometryFromJson(argv[2]);
#endif

          OGR_G_Segmentize(geometry, adfGeoTransform[1]);
          int nbGeom=OGR_G_GetPointCount(geometry);
          int k=0;
          double ppx=0,ppy=0;
          double value;
          char *buffer=NULL;
          int length=0;
          buffer=(char*)malloc(37*sizeof(char));
          sprintf(buffer,"{\"type\":\"LineString\",\"coordinates\":[");
          length+=strlen(buffer);
          for(k=0;k<nbGeom;k++){
            //OGRGeometryH point;
            double prx,pry,prz;
            OGR_G_GetPoint(geometry,k,&prx,&pry,&prz);
            float *pafScanline;
            pafScanline = (float *) CPLMalloc(sizeof(float));
            int px=(int)floor((prx-adfGeoTransform[0])/adfGeoTransform[1]);
            int py=(int)floor((pry-adfGeoTransform[3])/adfGeoTransform[5]);
            if(px!=ppx || py!=ppy){
              if(GDALRasterIO( hBand, GF_Read, px, py, 1, 1, 
                            pafScanline, 1, 1, GDT_Float32, 
                               0, 0 ) != CE_None){
                char *tmp;
                tmp=(char*) malloc(300*sizeof(char));
                sprintf(tmp,"GDALRasterIO failed for point (%d,%d)",px,py);
                setMapInMaps(conf,"lenv","message",_ss(tmp));
                CPLFree(pafScanline);
                free(tmp);
                return SERVICE_FAILED;
              }
              if(buffer!=NULL){
                int len=strlen(buffer);
                buffer=(char*)realloc(buffer,(len+50+1)*sizeof(char));
              }
              else
                buffer=(char*)malloc((51)*sizeof(char));
              char *tmpValue=(char *)malloc(50*sizeof(char));
              sprintf(tmpValue,"[%.6f,%.6f,%.6f]%c",prx,pry,pafScanline[0],(k+1==nbGeom?' ':','));
              strncpy(buffer+length,tmpValue,strlen(tmpValue));
              length+=strlen(tmpValue);
              buffer[length]=0;
              value=pafScanline[0];
              free(tmpValue);
              //Usefull if we can export 3D JSON string at the end
              //OGR_G_SetPoint(geometry,k,prx,pry,pafScanline[0]);            
            }
            else{
              if(buffer!=NULL)
                buffer=(char*)realloc(buffer,(strlen(buffer)+50+1)*sizeof(char));
              else
                buffer=(char*)malloc((51)*sizeof(char));
              char *tmpValue=(char *)malloc(50*sizeof(char));
              sprintf(tmpValue,"[%.6f,%.6f,%.6f]%c",prx,pry,value,(k+1==nbGeom?' ':','));
              strncpy(buffer+length,tmpValue,strlen(tmpValue));
              length+=strlen(tmpValue);
              buffer[length]=0;
              free(tmpValue);
              value=value;
            }
            CPLFree(pafScanline);
            ppx=px;
            ppy=py;
          }
          buffer=(char*)realloc(buffer,(strlen(buffer)+3)*sizeof(char));
          char *tmpValue=(char *)malloc(3*sizeof(char));
          sprintf(tmpValue,"]}");
          tmpValue[2]=0;
          strncpy(buffer+length,tmpValue,strlen(tmpValue));
          length+=strlen(tmpValue);
          buffer[length]=0;
#ifdef ZOO_SERVICE
          setMapInMaps(outputs,"Profile","value",buffer);
          setMapInMaps(outputs,"Profile","mimeType","text/plain");
#else
          fprintf(stderr,"%s\n",buffer);
#endif
          free(buffer);
          free(tmpValue);
          OGR_G_DestroyGeometry(geometry);
        }
    }
  else{
#ifdef ZOO_SERVICE
    setMapInMaps(conf,"lenv","message",_ss("Unable to load your raster file !"));
    return SERVICE_FAILED;
#else
    printf("Unable to load your raster file %s !\n",argv[1]);
#endif
  }
  GDALClose(hDataset);
  OGRCleanupAll();
  GDALDestroyDriverManager();
#ifdef ZOO_SERVICE
  return SERVICE_SUCCEEDED;
#endif
}

#ifdef ZOO_SERVICE
}
#endif