1 | /* **************************************************************************** |
---|
2 | * $Id$ |
---|
3 | * |
---|
4 | * Project: GdalExtractProfile |
---|
5 | * Purpose: Extract Profile from a Raster file for an Input Geometry (LINE) |
---|
6 | * Author: Gérald Fenoy, gerald.fenoy@geolabs.fr |
---|
7 | * |
---|
8 | * **************************************************************************** |
---|
9 | * Copyright (c) 2010, GeoLabs SARL |
---|
10 | * |
---|
11 | * Permission is hereby granted, free of charge, to any person obtaining a |
---|
12 | * copy of this software and associated documentation files (the "Software"), |
---|
13 | * to deal in the Software without restriction, including without limitation |
---|
14 | * the rights to use, copy, modify, merge, publish, distribute, sublicense, |
---|
15 | * and/or sell copies of the Software, and to permit persons to whom the |
---|
16 | * Software is furnished to do so, subject to the following conditions: |
---|
17 | * |
---|
18 | * The above copyright notice and this permission notice shall be included |
---|
19 | * in all copies or substantial portions of the Software. |
---|
20 | * |
---|
21 | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
---|
22 | * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
---|
23 | * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
---|
24 | * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
---|
25 | * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
---|
26 | * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
---|
27 | * DEALINGS IN THE SOFTWARE. |
---|
28 | ****************************************************************************/ |
---|
29 | |
---|
30 | #ifdef ZOO_SERVICE |
---|
31 | #include "service.h" |
---|
32 | #endif |
---|
33 | #include "gdal.h" |
---|
34 | #include "cpl_conv.h" |
---|
35 | #include "ogr_api.h" |
---|
36 | |
---|
37 | #ifdef ZOO_SERVICE |
---|
38 | extern "C" { |
---|
39 | #endif |
---|
40 | |
---|
41 | #ifdef ZOO_SERVICE |
---|
42 | int GdalExtractProfile(maps*& conf,maps*& inputs,maps*& outputs) |
---|
43 | #else |
---|
44 | int main(int argc,char** argv) |
---|
45 | #endif |
---|
46 | { |
---|
47 | char *pszFilename; |
---|
48 | #ifdef ZOO_SERVICE |
---|
49 | map* tmp=NULL; |
---|
50 | map* tmp1=NULL; |
---|
51 | tmp=getMapFromMaps(conf,"main","dataPath"); |
---|
52 | tmp1=getMapFromMaps(inputs,"RasterFile","value"); |
---|
53 | pszFilename=(char *)malloc((2+strlen(tmp->value)+strlen(tmp1->value))*sizeof(char)); |
---|
54 | sprintf(pszFilename,"%s/%s",tmp->value,tmp1->value); |
---|
55 | #else |
---|
56 | pszFilename=argv[1]; |
---|
57 | #endif |
---|
58 | GDALDatasetH hDataset; |
---|
59 | GDALAllRegister(); |
---|
60 | OGRRegisterAll(); |
---|
61 | |
---|
62 | hDataset = GDALOpen( pszFilename, GA_ReadOnly ); |
---|
63 | if( hDataset != NULL ) |
---|
64 | { |
---|
65 | GDALDriverH hDriver; |
---|
66 | double adfGeoTransform[6]; |
---|
67 | |
---|
68 | if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ) |
---|
69 | { |
---|
70 | |
---|
71 | |
---|
72 | GDALRasterBandH hBand; |
---|
73 | int nBlockXSize, nBlockYSize; |
---|
74 | int bGotMin, bGotMax; |
---|
75 | double adfMinMax[2]; |
---|
76 | |
---|
77 | hBand = GDALGetRasterBand( hDataset, 1 ); |
---|
78 | |
---|
79 | adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin ); |
---|
80 | adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax ); |
---|
81 | if( ! (bGotMin && bGotMax) ) |
---|
82 | GDALComputeRasterMinMax( hBand, TRUE, adfMinMax ); |
---|
83 | |
---|
84 | #ifdef ZOO_SERVICE |
---|
85 | tmp1=getMapFromMaps(inputs,"Geometry","value"); |
---|
86 | OGRGeometryH geometry=OGR_G_CreateGeometryFromJson(tmp1->value); |
---|
87 | #else |
---|
88 | OGRGeometryH geometry=OGR_G_CreateGeometryFromJson(argv[2]); |
---|
89 | #endif |
---|
90 | OGR_G_Segmentize(geometry, adfGeoTransform[1]); |
---|
91 | int nbGeom=OGR_G_GetPointCount(geometry); |
---|
92 | int k=0; |
---|
93 | double ppx=0,ppy=0; |
---|
94 | double value; |
---|
95 | char *buffer=NULL; |
---|
96 | int length=0; |
---|
97 | buffer=(char*)malloc(37*sizeof(char)); |
---|
98 | sprintf(buffer,"{\"type\":\"LineString\",\"coordinates\":["); |
---|
99 | length+=strlen(buffer); |
---|
100 | for(k=0;k<nbGeom;k++){ |
---|
101 | OGRGeometryH point; |
---|
102 | double prx,pry,prz; |
---|
103 | OGR_G_GetPoint(geometry,k,&prx,&pry,&prz); |
---|
104 | float *pafScanline; |
---|
105 | |
---|
106 | pafScanline = (float *) CPLMalloc(sizeof(float)); |
---|
107 | int px=(int)floor((prx-adfGeoTransform[0])/adfGeoTransform[1]); |
---|
108 | int py=(int)floor((pry-adfGeoTransform[3])/adfGeoTransform[5]); |
---|
109 | if(px!=ppx || py!=ppy){ |
---|
110 | if(GDALRasterIO( hBand, GF_Read, px, py, 1, 1, |
---|
111 | pafScanline, 1, 1, GDT_Float32, |
---|
112 | 0, 0 ) != CE_None){ |
---|
113 | char *tmp; |
---|
114 | tmp=(char*) malloc(300*sizeof(char)); |
---|
115 | sprintf(tmp,"GDALRasterIO failed for point (%d,%d)",px,py); |
---|
116 | setMapInMaps(conf,"lenv","message",_ss(tmp)); |
---|
117 | return SERVICE_FAILED; |
---|
118 | } |
---|
119 | if(buffer) |
---|
120 | buffer=(char*)realloc(buffer,(strlen(buffer)+50+1)*sizeof(char)); |
---|
121 | else |
---|
122 | buffer=(char*)malloc((51)*sizeof(char)); |
---|
123 | char *tmpValue=(char *)malloc(50*sizeof(char)); |
---|
124 | sprintf(tmpValue,"[%.6f,%.6f,%.6f]%c",prx,pry,pafScanline[0],(k+1==nbGeom?' ':',')); |
---|
125 | memcpy(buffer+length,tmpValue,strlen(tmpValue)); |
---|
126 | length+=strlen(tmpValue); |
---|
127 | value=pafScanline[0]; |
---|
128 | //Usefull if we can export 3D JSON string at the end |
---|
129 | //OGR_G_SetPoint(geometry,k,prx,pry,pafScanline[0]); |
---|
130 | } |
---|
131 | else{ |
---|
132 | if(buffer) |
---|
133 | buffer=(char*)realloc(buffer,(strlen(buffer)+50+1)*sizeof(char)); |
---|
134 | else |
---|
135 | buffer=(char*)malloc((51)*sizeof(char)); |
---|
136 | char *tmpValue=(char *)malloc(50*sizeof(char)); |
---|
137 | sprintf(tmpValue,"[%.6f,%.6f,%.6f]%c",prx,pry,value,(k+1==nbGeom?' ':',')); |
---|
138 | memcpy(buffer+length,tmpValue,strlen(tmpValue)); |
---|
139 | length+=strlen(tmpValue); |
---|
140 | value=value; |
---|
141 | } |
---|
142 | ppx=px; |
---|
143 | ppy=py; |
---|
144 | } |
---|
145 | buffer=(char*)realloc(buffer,(strlen(buffer)+3)*sizeof(char)); |
---|
146 | char *tmpValue=(char *)malloc(3*sizeof(char)); |
---|
147 | sprintf(tmpValue,"]}"); |
---|
148 | memcpy(buffer+length,tmpValue,strlen(tmpValue)); |
---|
149 | #ifdef ZOO_SERVICE |
---|
150 | setMapInMaps(outputs,"Profile","value",buffer); |
---|
151 | setMapInMaps(outputs,"Profile","mimeType","text/plain"); |
---|
152 | #else |
---|
153 | fprintf(stderr,"%s\n",buffer); |
---|
154 | #endif |
---|
155 | } |
---|
156 | } |
---|
157 | else{ |
---|
158 | #ifdef ZOO_SERVICE |
---|
159 | setMapInMaps(conf,"lenv","message",_ss("Unable to load your raster file !")); |
---|
160 | return SERVICE_FAILED; |
---|
161 | #else |
---|
162 | printf("Unable to load your raster file %s !\n",argv[1]); |
---|
163 | #endif |
---|
164 | } |
---|
165 | |
---|
166 | GDALClose(hDataset); |
---|
167 | #ifdef ZOO_SERVICE |
---|
168 | return SERVICE_SUCCEEDED; |
---|
169 | #endif |
---|
170 | } |
---|
171 | |
---|
172 | #ifdef ZOO_SERVICE |
---|
173 | } |
---|
174 | #endif |
---|