Skip to content

Commit

Permalink
Removed old lon offset code
Browse files Browse the repository at this point in the history
  • Loading branch information
ScottMcMichael committed May 3, 2016
1 parent 69fd948 commit 6f0f1b7
Show file tree
Hide file tree
Showing 5 changed files with 7 additions and 63 deletions.
57 changes: 5 additions & 52 deletions src/vw/Cartography/GeoTransform.cc
Expand Up @@ -36,43 +36,7 @@ namespace cartography {
// Constructor
GeoTransform::GeoTransform(GeoReference const& src_georef, GeoReference const& dst_georef,
BBox2i const& src_bbox, BBox2i const& dst_bbox) :
m_src_georef(src_georef), m_dst_georef(dst_georef), m_offset(Vector2(0, 0)) {

// Deal with the fact that longitudes could differ by 360 degrees

// between src and dst.
if (src_bbox.empty() || dst_bbox.empty()) {
// If we don't know the image areas, simply use the 0, 0 corner
try {
Vector2 src_origin = src_georef.pixel_to_lonlat(Vector2(0, 0));
Vector2 dst_origin = dst_georef.pixel_to_lonlat(Vector2(0, 0));
m_offset = Vector2( 360.0*round( (dst_origin[0] - src_origin[0])/360.0 ), 0.0 );
}catch ( const std::exception & e ) {
m_offset = Vector2(0, 0);
}
} else{

// Try to offset by 360 degrees until the lon-lat boxes are most compatible,
// measured by how close the centers of the boxes are.

BBox2 src_lonlat_box = src_georef.pixel_to_lonlat_bbox(src_bbox);
BBox2 dst_lonlat_box = dst_georef.pixel_to_lonlat_bbox(dst_bbox);
std::vector<double> shift, dist_vec;
for (int val = -360; val <= 360; val+= 360) {
shift.push_back(val);
BBox2 shifted_src = src_lonlat_box + Vector2(val, 0);

// Distance between centers of boxes. The smaller it is, the better.
double dist = norm_2((shifted_src.min() + shifted_src.max())/2.0
- (dst_lonlat_box.min() + dst_lonlat_box.max())/2.0);
dist_vec.push_back(dist);
}
int min_index = std::distance(dist_vec.begin(), min_element(dist_vec.begin(), dist_vec.end()));
m_offset = Vector2(shift[min_index], 0);
}

//std::cout << "Geotransform: m_offset = " << m_offset << std::endl;
m_offset = Vector2(0, 0); // Disable offset code!
m_src_georef(src_georef), m_dst_georef(dst_georef) {

const std::string src_datum = m_src_georef.datum().proj4_str();
const std::string dst_datum = m_dst_georef.datum().proj4_str();
Expand All @@ -85,14 +49,6 @@ namespace cartography {
else
m_skip_map_projection = false;

// Bugfix: If one of the projections is lon-lat, and the offset is 360 degrees,
// don't skip map-projection, as have to correct for the offset.
bool has_lonat = (m_src_georef.proj4_str().find("+proj=longlat") != std::string::npos ||
m_dst_georef.proj4_str().find("+proj=longlat") != std::string::npos );
if (has_lonat && m_offset != Vector2()) {
m_skip_map_projection = false;
}

// This optimizes the case where the two datums are the same,
// and thus we don't need to call proj to convert between them
// as we transform.
Expand All @@ -119,14 +75,11 @@ namespace cartography {
set_tolerance( 0.1 );
}

void GeoTransform::set_offset(Vector2 const& offset){
m_offset = offset;
}

Vector2 GeoTransform::reverse(Vector2 const& v) const {
if (m_skip_map_projection)
return m_src_georef.point_to_pixel(m_dst_georef.pixel_to_point(v));
Vector2 src_lonlat = m_dst_georef.pixel_to_lonlat(v) - m_offset;
Vector2 src_lonlat = m_dst_georef.pixel_to_lonlat(v);
if(m_skip_datum_conversion)
return m_src_georef.lonlat_to_pixel(src_lonlat);
src_lonlat = lonlat_to_lonlat(src_lonlat, false);
Expand All @@ -137,7 +90,7 @@ namespace cartography {
Vector2 GeoTransform::forward(Vector2 const& v) const {
if (m_skip_map_projection)
return m_dst_georef.point_to_pixel(m_src_georef.pixel_to_point(v));
Vector2 dst_lonlat = m_src_georef.pixel_to_lonlat(v) + m_offset;
Vector2 dst_lonlat = m_src_georef.pixel_to_lonlat(v);
if(m_skip_datum_conversion)
return m_dst_georef.lonlat_to_pixel(dst_lonlat);
dst_lonlat = lonlat_to_lonlat(dst_lonlat, true);
Expand Down Expand Up @@ -192,7 +145,7 @@ namespace cartography {
Vector2 GeoTransform::pixel_to_point( Vector2 const& v ) const {
if (m_skip_map_projection)
return m_src_georef.pixel_to_point(v);
Vector2 src_lonlat = m_src_georef.pixel_to_lonlat(v) + m_offset;
Vector2 src_lonlat = m_src_georef.pixel_to_lonlat(v);
if(m_skip_datum_conversion)
return m_dst_georef.lonlat_to_point(src_lonlat);
Vector2 dst_lonlat = lonlat_to_lonlat(src_lonlat, true);
Expand All @@ -202,7 +155,7 @@ namespace cartography {
Vector2 GeoTransform::point_to_pixel( Vector2 const& v ) const {
if (m_skip_map_projection)
return m_src_georef.point_to_pixel(v);
Vector2 src_lonlat = m_src_georef.point_to_lonlat(v) + m_offset;
Vector2 src_lonlat = m_src_georef.point_to_lonlat(v);
if(m_skip_datum_conversion)
return m_dst_georef.lonlat_to_pixel(src_lonlat);
Vector2 dst_lonlat = lonlat_to_lonlat(src_lonlat, true);
Expand Down
5 changes: 0 additions & 5 deletions src/vw/Cartography/GeoTransform.h
Expand Up @@ -44,7 +44,6 @@ namespace cartography {
m_dst_datum_proj;
bool m_skip_map_projection;
bool m_skip_datum_conversion;
Vector2 m_offset; // TODO delete

public:
/// Normal constructor
Expand Down Expand Up @@ -102,10 +101,6 @@ namespace cartography {
/// Convert a point bounding box in the source to a pixel bounding box in the destination.
BBox2 point_to_pixel_bbox( BBox2 const& point_bbox ) const;


/// Sometimes the offset is not computed correctly, for example, when one
/// of the two georeferenced images encompasses the whole globe.
void set_offset(Vector2 const& offset);

}; // End class GeoTransform

Expand Down
2 changes: 1 addition & 1 deletion src/vw/Cartography/tests/TestGeoReferenceUtils.cxx
Expand Up @@ -97,7 +97,7 @@ TEST( Common, gdal_write_checks ) {
double nodata = -1000;
bool has_nodata = true, has_georef = true;
TerminalProgressCallback tpc("vw", "");
BaseOptions opt;
GdalWriteOptions opt;

block_write_gdal_image("img1.tif", dem, has_georef, georef, has_nodata, nodata, opt, tpc);
block_write_gdal_image("img2.tif", dem, has_georef, georef, has_nodata, nodata, opt, tpc);
Expand Down
2 changes: 1 addition & 1 deletion src/vw/Cartography/tests/TestGeoTransform.cxx
Expand Up @@ -147,7 +147,7 @@ TEST(GeoTransform, RefToRef) {
//std::cout << "pixel1 = " << pixel1 << std::endl;
//std::cout << "point1 = " << point1 << std::endl;

Vector2 lonlat = georef1.point_to_lonlat(point1);
//Vector2 lonlat = georef1.point_to_lonlat(point1);
//Vector2 point2 = georef2.lonlat_to_point(lonlat);
//Vector2 pixel2 = georef2.point_to_pixel(point2);
//std::cout << "lonlat = " << lonlat << std::endl;
Expand Down
4 changes: 0 additions & 4 deletions src/vw/tools/image2qtree.h
Expand Up @@ -274,10 +274,6 @@ void do_mosaic(const Options& opt, const vw::ProgressCallback *progress) {
boost::shared_ptr<DiskImageResource> file( DiskImageResource::open(filename) );
GeoTransform geotx( input_georef, output_georef );

// Even though the output georeference starts at -180, and input georef may start
// close to 180, here we don't want to correct for this discrepancy.
geotx.set_offset(Vector2(0, 0));

ImageViewRef<PixelT> source = DiskImageView<PixelT>( file );

// Handle nodata values/mask
Expand Down

0 comments on commit 6f0f1b7

Please sign in to comment.