Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

WIP

  • Loading branch information...
commit 56bded5f4b6e356f5757e2f946a6558304e3aa37 1 parent b50ff49
Zack Moratto authored September 02, 2012
54  src/vw/Cartography/tests/TestGeoTransform.cxx
@@ -23,6 +23,10 @@
23 23
 #include <vw/Cartography/GeoReference.h>
24 24
 #include <vw/Cartography/GeoTransform.h>
25 25
 
  26
+#if defined(VW_HAVE_PKG_GDAL) && VW_HAVE_PKG_GDAL==1
  27
+#include "ogr_spatialref.h"
  28
+#endif
  29
+
26 30
 using namespace vw;
27 31
 using namespace vw::cartography;
28 32
 using namespace vw::test;
@@ -108,3 +112,53 @@ TEST( GeoTransform, StereographicSingularity ) {
108 112
   EXPECT_NO_THROW( output = geotx.forward_bbox(input) );
109 113
   EXPECT_NEAR( 0, output.min()[1], 2 );
110 114
 }
  115
+
  116
+void load_proj4( std::string const& proj4,
  117
+                 GeoReference& georef ) {
  118
+  OGRSpatialReference gdal_spatial_ref;
  119
+  if (gdal_spatial_ref.SetFromUserInput( proj4.c_str() ))
  120
+    vw_throw( ArgumentErr() << "Failed to parse: \"" << proj4 << "\"." );
  121
+  char *wkt = NULL;
  122
+  gdal_spatial_ref.exportToWkt( &wkt );
  123
+  std::string wkt_string(wkt);
  124
+  delete[] wkt;
  125
+
  126
+  georef.set_wkt( wkt_string );
  127
+}
  128
+
  129
+TEST( GeoTransform, DatumTransition ) {
  130
+  GeoReference georef1, georef2, georef3;
  131
+  load_proj4("+proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=609601.2192024384 +y_0=0 +ellps=clrk66 +datum=NAD27 +to_meter=0.34 +no_defs", georef1);
  132
+  georef1.set_transform(math::identity_matrix<3>());
  133
+  load_proj4("+proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",georef2);
  134
+  Matrix3x3 transform = math::identity_matrix<3>();
  135
+  transform(0,0) = 0.5;
  136
+  transform(1,1) = -0.5;
  137
+  georef2.set_transform(math::identity_matrix<3>());
  138
+  load_proj4("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ",georef3);
  139
+  transform(0,2) = 10000;
  140
+  transform(1,2) = -1000000;
  141
+  georef3.set_transform(transform);
  142
+
  143
+  GeoReference georef2_ll;
  144
+  load_proj4("+proj=longlat +units=m",georef2_ll);
  145
+  georef2_ll.set_transform(math::identity_matrix<3>());
  146
+
  147
+  // Comparing 1->2 to 1->3->2. They should produce the same results if everything worked
  148
+  GeoTransform geotx12( georef1, georef2 );
  149
+  GeoTransform geotx13( georef1, georef3 );
  150
+  GeoTransform geotx32( georef3, georef2 );
  151
+
  152
+  for ( size_t x = 0; x < 10000000; x += 100000 ) {
  153
+    for ( size_t y = 0; y < 10000000; y+= 100000 ) {
  154
+      // Verify GeoTransforms are consistent
  155
+      EXPECT_VECTOR_NEAR( geotx12.forward( Vector2(x,y) ),
  156
+                          geotx32.forward(geotx13.forward( Vector2(x,y) )), 1e-6 );
  157
+
  158
+      // Verify that we're consistent with GeoRefs
  159
+      Vector2 lonlat1 = georef1.pixel_to_lonlat( Vector2(x,y) ); // it appears like lonlat is in geocentric coordinates.
  160
+      EXPECT_VECTOR_NEAR( geotx12.forward( Vector2(x,y) ),
  161
+                          georef2.lonlat_to_pixel(lonlat1), 1e-6 );
  162
+    }
  163
+  }
  164
+}

0 notes on commit 56bded5

Please sign in to comment.
Something went wrong with that request. Please try again.