← Back to branch summary

~ubuntu-branches/ubuntu/hardy/qgis/hardy

« back to all changes in this revision

Viewing changes to src/qgscoordinatetransform.cpp

  • Committer: Bazaar Package Importer
  • Author(s): William Grant
  • Date: 2007-05-06 13:42:32 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20070506134232-pyli6t388w5asd8x
Tags: 0.8.0-3ubuntu1
* Merge from Debian unstable. Remaining Ubuntu changes:
  - debian/rules, debian/qgis.install, debian/qgis.dirs debian/qgis.desktop:
    Add and install .desktop.
* debian/qgis.desktop: Remove Applications category; it's not real.
* Modify Maintainer value to match Debian-Maintainer-Field Spec

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/qgscoordinatetransform.cpp
1
 
/***************************************************************************
2
 
               QgsCoordinateTransform.cpp  - Coordinate Transforms
3
 
                             -------------------
4
 
    begin                : Dec 2004
5
 
    copyright            : (C) 2004 Tim Sutton
6
 
    email                : tim at linfiniti.com
7
 
 ***************************************************************************/
8
 
 
9
 
/***************************************************************************
10
 
 *                                                                         *
11
 
 *   This program is free software; you can redistribute it and/or modify  *
12
 
 *   it under the terms of the GNU General Public License as published by  *
13
 
 *   the Free Software Foundation; either version 2 of the License, or     *
14
 
 *   (at your option) any later version.                                   *
15
 
 *                                                                         *
16
 
 ***************************************************************************/
17
 
/* $Id: qgscoordinatetransform.cpp,v 1.42.2.16 2005/09/21 22:51:25 timlinux Exp $ */
18
 
#include <cassert>
19
 
#include "qgscoordinatetransform.h"
20
 
 
21
 
//qt includes
22
 
#include <qdom.h>
23
 
 
24
 
QgsCoordinateTransform::QgsCoordinateTransform( ) : QObject()
25
 
 
26
 
{
27
 
}
28
 
 
29
 
QgsCoordinateTransform::QgsCoordinateTransform(const QgsSpatialRefSys& source, 
30
 
                                               const QgsSpatialRefSys& dest)
31
 
{
32
 
  mSourceSRS = source;
33
 
  mDestSRS = dest;
34
 
  initialise();
35
 
}
36
 
 
37
 
 
38
 
QgsCoordinateTransform::QgsCoordinateTransform( QString theSourceSRS, QString theDestSRS ) : QObject()
39
 
 
40
 
{
41
 
  mSourceSRS.createFromWkt(theSourceSRS);
42
 
  mDestSRS.createFromWkt(theDestSRS);
43
 
  // initialize the coordinate system data structures
44
 
  //XXX Who spells initialize initialise?
45
 
  //XXX A: Its the queen's english....
46
 
  //XXX  : Long live the queen! Lets get on with the initialisation...
47
 
  initialise();
48
 
}
49
 
 
50
 
QgsCoordinateTransform::QgsCoordinateTransform(long theSourceSrid,
51
 
    QString theDestWKT,
52
 
    QgsSpatialRefSys::SRS_TYPE theSourceSRSType): QObject()
53
 
{
54
 
 
55
 
  mSourceSRS.createFromId(theSourceSrid, theSourceSRSType);
56
 
  mDestSRS.createFromWkt(theDestWKT);
57
 
  // initialize the coordinate system data structures
58
 
  //XXX Who spells initialize initialise?
59
 
  //XXX A: Its the queen's english....
60
 
  //XXX  : Long live the queen! Lets get on with the initialisation...
61
 
  initialise();
62
 
}
63
 
 
64
 
QgsCoordinateTransform::~QgsCoordinateTransform()
65
 
{
66
 
  // free the proj objects
67
 
  if (mSourceProjection) 
68
 
  {
69
 
    pj_free(mSourceProjection);
70
 
  }
71
 
  if (mDestinationProjection)
72
 
  {
73
 
    pj_free(mDestinationProjection);
74
 
  }
75
 
}
76
 
 
77
 
void QgsCoordinateTransform::setSourceSRS(const QgsSpatialRefSys& theSRS)
78
 
{
79
 
  mSourceSRS = theSRS;
80
 
  initialise();
81
 
}
82
 
void QgsCoordinateTransform::setDestSRS(const QgsSpatialRefSys& theSRS)
83
 
{
84
 
#ifdef QGISDEBUG
85
 
  std::cout << "QgsCoordinateTransform::setDestSRS called" << std::endl;
86
 
#endif
87
 
  mDestSRS = theSRS;
88
 
  initialise();
89
 
}
90
 
 
91
 
 
92
 
void QgsCoordinateTransform::setDestSRSID (long theSRSID)
93
 
{
94
 
  //!todo Add some logic here to determine if the srsid is a system or user one
95
 
#ifdef QGISDEBUG
96
 
  std::cout << "QgsCoordinateTransform::setDestSRSID slot called" << std::endl;
97
 
#endif
98
 
  mDestSRS.createFromSrsId(theSRSID);
99
 
  initialise();
100
 
}
101
 
 
102
 
// XXX This whole function is full of multiple return statements!!!
103
 
void QgsCoordinateTransform::initialise()
104
 
{
105
 
 
106
 
  mInitialisedFlag=false; //guilty until proven innocent...
107
 
  mSourceProjection = NULL;
108
 
  mDestinationProjection = NULL;
109
 
 
110
 
  // XXX Warning - multiple return paths in this block!!
111
 
  if (!mSourceSRS.isValid())
112
 
  {
113
 
    //mSourceSRS = defaultWkt;
114
 
    // Pass through with no projection since we have no idea what the layer
115
 
    // coordinates are and projecting them may not be appropriate
116
 
    mShortCircuit = true;
117
 
    return;
118
 
  }
119
 
 
120
 
  if (!mDestSRS.isValid())
121
 
  {
122
 
    //No destination projection is set so we set the default output projection to
123
 
    //be the same as input proj. This only happens on the first layer loaded
124
 
    //whatever that may be...
125
 
    mDestSRS.createFromProj4(mSourceSRS.proj4String());
126
 
  }
127
 
 
128
 
  //XXX todo overload == operator for QgsSpatialRefSys
129
 
  //at the moment srs.parameters contains the whole proj def...soon it wont...
130
 
  //if (mSourceSRS->proj4String() == mDestSRS->proj4String())
131
 
  if (mSourceSRS == mDestSRS)
132
 
  {
133
 
    // If the source and destination projection are the same, set the short
134
 
    // circuit flag (no transform takes place)
135
 
    mShortCircuit=true;
136
 
    return;
137
 
  }
138
 
  else
139
 
  {
140
 
    // Transform must take place
141
 
    mShortCircuit=false;
142
 
  }
143
 
 
144
 
  // init the projections (destination and source)
145
 
  mDestinationProjection = pj_init_plus(mDestSRS.proj4String().local8Bit());
146
 
  mSourceProjection = pj_init_plus(mSourceSRS.proj4String().local8Bit());
147
 
 
148
 
  mInitialisedFlag = true;
149
 
  if ( mDestinationProjection == NULL )
150
 
  {
151
 
    mInitialisedFlag = false;
152
 
  }
153
 
  if ( mSourceProjection == NULL )
154
 
  {
155
 
    mInitialisedFlag = false;
156
 
  }
157
 
 
158
 
  if (mInitialisedFlag)
159
 
  {
160
 
 
161
 
    std::cout << "------------------------------------------------------------"<< std::endl;
162
 
    std::cout << "QgsCoordinateTransform::initialise()" << std::endl;
163
 
    std::cout << "The OGR Coordinate transformation for this layer was set to" << std::endl;
164
 
    // note overloaded << operator on qgsspatialrefsys cant be used on pointers -
165
 
    // so we dereference them like this (*mSourceSRS) (Thanks Lars for pointing that out)
166
 
    std::cout << "INPUT: " << std::endl << mSourceSRS << std::endl;
167
 
    std::cout << "OUTPUT: " << std::endl << mDestSRS  << std::endl;
168
 
    std::cout << "------------------------------------------------------------" << std::endl;
169
 
  }
170
 
  else
171
 
  {
172
 
    std::cout << "------------------------------------------------------------"<< std::endl;
173
 
    std::cout << "QgsCoordinateTransform::initialise()" << std::endl;
174
 
    std::cout << "The OGR Coordinate transformation FAILED TO INITIALISE!" << std::endl;
175
 
    std::cout << "------------------------------------------------------------"<< std::endl;
176
 
  }
177
 
}
178
 
 
179
 
//
180
 
//
181
 
// TRANSFORMERS BELOW THIS POINT .........
182
 
//
183
 
//
184
 
//
185
 
 
186
 
 
187
 
QgsPoint QgsCoordinateTransform::transform(const QgsPoint thePoint,TransformDirection direction)
188
 
{
189
 
  if (mShortCircuit || !mInitialisedFlag) return thePoint;
190
 
  // transform x
191
 
  double x = thePoint.x();
192
 
  double y = thePoint.y();
193
 
  double z = 0.0;
194
 
  try
195
 
  { 
196
 
    transformCoords(1, &x, &y, &z, direction );
197
 
  }
198
 
  catch(QgsCsException &cse)
199
 
  {
200
 
    // rethrow the exception
201
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
202
 
    throw cse;
203
 
  }
204
 
 
205
 
#ifdef QGISDEBUG
206
 
  //std::cout << "Point projection...X : " << thePoint.x() << "-->" << x << ", Y: " << thePoint.y() << " -->" << y << std::endl;
207
 
#endif
208
 
  return QgsPoint(x, y);
209
 
}
210
 
 
211
 
 
212
 
QgsPoint QgsCoordinateTransform::transform(const double theX, const double theY=0,TransformDirection direction)
213
 
{
214
 
  try
215
 
  {
216
 
    return transform(QgsPoint(theX, theY), direction);
217
 
  }
218
 
  catch(QgsCsException &cse)
219
 
  {
220
 
    // rethrow the exception
221
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
222
 
    throw cse;
223
 
  }
224
 
}
225
 
 
226
 
QgsRect QgsCoordinateTransform::transform(const QgsRect theRect,TransformDirection direction)
227
 
{
228
 
  if (mShortCircuit || !mInitialisedFlag) return theRect;
229
 
  // transform x
230
 
  double x1 = theRect.xMin();
231
 
  double y1 = theRect.yMin();
232
 
  double x2 = theRect.xMax();
233
 
  double y2 = theRect.yMax();
234
 
 
235
 
#ifdef QGISDEBUG
236
 
  std::cout << this;
237
 
#endif
238
 
  // Number of points to reproject------+
239
 
  //                                    |
240
 
  //                                    V
241
 
  try
242
 
  {
243
 
    double z = 0.0;
244
 
    transformCoords(1, &x1, &y1, &z, direction);
245
 
    transformCoords(1, &x2, &y2, &z, direction);
246
 
  }
247
 
  catch(QgsCsException &cse)
248
 
  {
249
 
    // rethrow the exception
250
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
251
 
    throw cse;
252
 
  }
253
 
 
254
 
#ifdef QGISDEBUG
255
 
  std::cout << "Rect projection..."
256
 
  << "Xmin : "
257
 
  << theRect.xMin()
258
 
  << "-->" << x1
259
 
  << ", Ymin: "
260
 
  << theRect.yMin()
261
 
  << " -->" << y1
262
 
  << "Xmax : "
263
 
  << theRect.xMax()
264
 
  << "-->" << x2
265
 
  << ", Ymax: "
266
 
  << theRect.yMax()
267
 
  << " -->" << y2
268
 
  << std::endl;
269
 
#endif
270
 
  return QgsRect(x1, y1, x2 , y2);
271
 
}
272
 
 
273
 
void QgsCoordinateTransform::transformInPlace(double& x, double& y, double& z,
274
 
    TransformDirection direction)
275
 
{
276
 
  if (mShortCircuit || !mInitialisedFlag)
277
 
    return;
278
 
#ifdef QGISDEBUG
279
 
  //std::cout << "Using transform in place " << __FILE__ << " " << __LINE__ << std::endl;
280
 
#endif
281
 
  // transform x
282
 
  try
283
 
  {
284
 
    transformCoords(1, &x, &y, &z, direction );
285
 
  }
286
 
  catch(QgsCsException &cse)
287
 
  {
288
 
    // rethrow the exception
289
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
290
 
    throw cse;
291
 
  }
292
 
}
293
 
 
294
 
void QgsCoordinateTransform::transformInPlace(std::vector<double>& x,
295
 
    std::vector<double>& y, std::vector<double>& z,
296
 
    TransformDirection direction)
297
 
{
298
 
  if (mShortCircuit || !mInitialisedFlag)
299
 
    return;
300
 
 
301
 
  assert(x.size() == y.size());
302
 
 
303
 
  // Apparently, if one has a std::vector, it is valid to use the
304
 
  // address of the first element in the vector as a pointer to an
305
 
  // array of the vectors data, and hence easily interface with code
306
 
  // that wants C-style arrays.
307
 
 
308
 
  try
309
 
  {
310
 
    transformCoords(x.size(), &x[0], &y[0], &z[0], direction);
311
 
  }
312
 
  catch(QgsCsException &cse)
313
 
  {
314
 
    // rethrow the exception
315
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
316
 
    throw cse;
317
 
  }
318
 
}
319
 
 
320
 
 
321
 
QgsRect QgsCoordinateTransform::transformBoundingBox(const QgsRect rect, TransformDirection direction)
322
 
{
323
 
  // Calculate the bounding box of a QgsRect in the source SRS
324
 
  // when projected to the destination SRS (or the inverse).
325
 
  // This is done by looking at a number of points spread evenly
326
 
  // across the rectangle
327
 
 
328
 
  if (mShortCircuit || !mInitialisedFlag)
329
 
    return rect;
330
 
 
331
 
  static const int numP = 8;
332
 
 
333
 
  QgsRect bb_rect;
334
 
  bb_rect.setMinimal();
335
 
 
336
 
  // We're interfacing with C-style vectors in the
337
 
  // end, so let's do C-style vectors here too.
338
 
  
339
 
  double x[numP * numP];
340
 
  double y[numP * numP];
341
 
  double z[numP * numP];
342
 
#ifdef QGISDEBUG   
343
 
  std::cout << "Entering transformBoundingBox..." << std::endl;
344
 
#endif 
345
 
  // Populate the vectors
346
 
 
347
 
  double dx = rect.width()  / (double)(numP - 1);
348
 
  double dy = rect.height() / (double)(numP - 1);
349
 
 
350
 
  double pointY = rect.yMin();
351
 
 
352
 
  for (int i = 0; i < numP ; i++)
353
 
  {
354
 
 
355
 
    // Start at right edge
356
 
    double pointX = rect.xMin();
357
 
 
358
 
    for (int j = 0; j < numP; j++)
359
 
    {
360
 
      x[(i*numP) + j] = pointX;
361
 
      y[(i*numP) + j] = pointY;
362
 
      // and the height...
363
 
      z[(i*numP) + j] = 0.0;
364
 
//      std::cout << "BBox coord: (" << x[(i*numP) + j] << ", " << y[(i*numP) + j]
365
 
//                << ")" << std::endl;
366
 
      pointX += dx; 
367
 
    }
368
 
    pointY += dy;
369
 
  }
370
 
 
371
 
  // Do transformation. Any exception generated must
372
 
  // be handled in above layers.
373
 
  try
374
 
  {
375
 
    transformCoords(numP * numP, x, y, z, direction);
376
 
  }
377
 
  catch(QgsCsException &cse)
378
 
  {
379
 
  // rethrow the exception
380
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
381
 
    throw cse;
382
 
  }
383
 
 
384
 
  // Calculate the bounding box and use that for the extent
385
 
        
386
 
  for (int i = 0; i < numP * numP; i++)
387
 
  {
388
 
    bb_rect.combineExtentWith(x[i], y[i]);
389
 
  }
390
 
#ifdef QGISDEBUG 
391
 
  std::cout << "Projected extent: " << (bb_rect.stringRep()).local8Bit() << std::endl;
392
 
#endif 
393
 
  return bb_rect;
394
 
}
395
 
 
396
 
void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z,TransformDirection direction)
397
 
{
398
 
  assert(mSourceSRS.isValid());
399
 
  assert(mDestSRS.isValid());
400
 
#ifdef QGISDEBUG
401
 
  //double xorg = x;
402
 
  //double yorg = y;
403
 
  //std::cout << "[[[[[[Number of points to transform: " << numPoints << "]]]]]]" << std::endl;
404
 
#endif
405
 
  // use proj4 to do the transform
406
 
  QString dir;
407
 
  // if the source/destination projection is lat/long, convert the points to radians
408
 
  // prior to transforming
409
 
  if((pj_is_latlong(mDestinationProjection) && (direction == INVERSE))
410
 
      || (pj_is_latlong(mSourceProjection) && (direction == FORWARD)))
411
 
  {
412
 
    for (int i = 0; i < numPoints; ++i)
413
 
    {
414
 
      x[i] *= DEG_TO_RAD;
415
 
      y[i] *= DEG_TO_RAD;
416
 
      z[i] *= DEG_TO_RAD;
417
 
    }
418
 
 
419
 
  }
420
 
  int projResult;
421
 
  if(direction == INVERSE)
422
 
  {
423
 
    /*
424
 
    std::cout << "!!!! INVERSE PROJ4 TRANSFORM !!!!" << std::endl; 
425
 
    std::cout << "     numPoint: " << numPoints << std::endl; 
426
 
    std::cout << "     x       : " << x << std::endl; 
427
 
    std::cout << "     y       : " << y << std::endl; 
428
 
    */
429
 
    projResult = pj_transform(mDestinationProjection, mSourceProjection , numPoints, 0, x, y, z);
430
 
    dir = "inverse";
431
 
  }
432
 
  else
433
 
  {
434
 
    /*
435
 
    std::cout << "!!!! FORWARD PROJ4 TRANSFORM !!!!" << std::endl; 
436
 
    std::cout << "     numPoint: " << numPoints << std::endl; 
437
 
    std::cout << "     x       : " << x << std::endl; 
438
 
    std::cout << "     y       : " << y << std::endl; 
439
 
    std::cout << "     z       : " << z << std::endl; 
440
 
    */
441
 
    assert(mSourceProjection != 0);
442
 
    assert(mDestinationProjection !=0);
443
 
    projResult = pj_transform(mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z);
444
 
    dir = "forward";
445
 
  }
446
 
 
447
 
  if (projResult != 0)
448
 
  {
449
 
    //something bad happened....
450
 
    QString msg;
451
 
    QTextOStream pjErr(&msg);
452
 
 
453
 
    pjErr << tr("Failed") << " " << dir << " " << tr("transform of") << '\n';
454
 
    for (int i = 0; i < numPoints; ++i)
455
 
    {
456
 
      if(direction == FORWARD)
457
 
      {
458
 
        pjErr << "(" << x[i] << ", " << y[i] << ")\n";
459
 
      }
460
 
      else
461
 
      {
462
 
        pjErr << "(" << x[i] * RAD_TO_DEG << ", " << y[i] * RAD_TO_DEG << ")\n";
463
 
      }
464
 
    }
465
 
 
466
 
    pjErr << tr("with error: ") << pj_strerrno(projResult) << '\n';
467
 
#ifdef QGISDEBUG
468
 
  std::cout << "Projection failed emitting invalid transform signal: \n" << msg.local8Bit() << std::endl;
469
 
#endif
470
 
    emit invalidTransformInput();
471
 
    std::cout << "Throwing exception " << __FILE__ << __LINE__ << std::endl; 
472
 
    throw  QgsCsException(msg);
473
 
  }
474
 
  // if the result is lat/long, convert the results from radians back
475
 
  // to degrees
476
 
  if((pj_is_latlong(mDestinationProjection) && (direction == FORWARD))
477
 
      || (pj_is_latlong(mSourceProjection) && (direction == INVERSE)))
478
 
  {
479
 
    for (int i = 0; i < numPoints; ++i)
480
 
    {
481
 
      x[i] *= RAD_TO_DEG;
482
 
      y[i] *= RAD_TO_DEG;
483
 
      z[i] *= RAD_TO_DEG;
484
 
    }
485
 
  }
486
 
#ifdef QGISDEBUG
487
 
  // std::cout << "[[[[[[ Projected " << xorg << ", " << yorg << " to "  << x << ", " << y << " ]]]]]]"<< std::endl;
488
 
#endif
489
 
}
490
 
 
491
 
bool QgsCoordinateTransform::readXML( QDomNode & theNode )
492
 
{
493
 
#ifdef QGISDEBUG
494
 
  std::cout << "Reading Coordinate Transform from xml ------------------------!" << std::endl;
495
 
#endif
496
 
  QDomNode mySrcNodeParent = theNode.namedItem("sourcesrs");
497
 
  QDomNode mySrcNode = mySrcNodeParent.namedItem("spatialrefsys");
498
 
  mSourceSRS.readXML(mySrcNode);
499
 
  QDomNode myDestNodeParent = theNode.namedItem("destinationsrs");
500
 
  QDomNode myDestNode = myDestNodeParent.namedItem("spatialrefsys");
501
 
  mDestSRS.readXML(myDestNode);
502
 
  initialise();
503
 
#ifdef WIN32
504
 
  return true;
505
 
#endif
506
 
}
507
 
 
508
 
bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
509
 
{
510
 
  
511
 
  QDomElement myNodeElement = theNode.toElement();
512
 
  QDomElement myTransformElement  = theDoc.createElement( "coordinatetransform" );
513
 
  
514
 
  QDomElement mySourceElement  = theDoc.createElement( "sourcesrs" );
515
 
  mSourceSRS.writeXML(mySourceElement, theDoc);
516
 
  myTransformElement.appendChild(mySourceElement);
517
 
  
518
 
  QDomElement myDestElement  = theDoc.createElement( "destinationsrs" );
519
 
  mDestSRS.writeXML(myDestElement, theDoc);
520
 
  myTransformElement.appendChild(myDestElement);
521
 
  
522
 
  myNodeElement.appendChild(myTransformElement);
523
 
#ifdef WIN32
524
 
  return true;
525
 
#endif
526
 
}
527