Opened 19 years ago

Last modified 19 years ago

#642 closed defect (fixed)

latlong to ortho projection fails

Reported by: tylermitchell@… Owned by: warmerdam
Priority: high Milestone:
Component: GDAL_Raster Version: unspecified
Severity: major Keywords:
Cc:

Description

Link to image in URL. 
----------- 
Re: [Gdal-dev] gdalwarp failing on proj=ortho 
 From: Frank Warmerdam <warmerdam@pobox.com> 
 To: Tyler Mitchell <tylermitchell@shaw.ca> 
 CC: gdal-dev@xserve.flids.com 
  
Tyler Mitchell wrote: 
> Sorry to bring up orthographic projections again, I know you looked into 
these  
> earlier this year.  I really want to show this one off though!  I'm trying 
to  
> convert a lat/long image into ortho - but I think I'm missing something. 
>  
> Here is my command line: 
>  
> gdalwarp -s_srs '+proj=latlong'  
>       -t_srs '+proj=ortho +ellps=WGS84 +lat_0=60 +lon_0=-95'  
>       globe.tif globe_ortho.tif 
>  
> The output error I get is:     
>  
> ERROR 1: tolerance condition error 
> ERROR 1: GDALWarperOperation::ComputeSourceWindow() failed because 
> the pfnTransformer failed. 
>  
> I'm running proj-4.4.7 and gdal from cvs a couple weeks ago.  Any ideas why  
> this doesn't work? 
 
Tyler, 
 
The tolerence condition error is coming out of PROJ.4. I don't know exactly 
what it means, but I presume it indicates some sort of mathematical error 
in transforming "over the horizon" or off the earth. 
 
If you can make globe.tif available to me I can see why gdalwarp isn't more 
adaptive when only subregions of the target area support reprojection.  I 
did a bunch of work on global scale transformations in MapServer but I don't 
think the same capabilities were implemented in GDAL itself.  In fact, I 
would appreciate it if you could submit this problem via the GDAL bugzilla. 
 
I would add that you should specify an ellipe or datum with all projections 
that don't implicitly contain one.  For instance you should have 
'+proj=latlong +ellps=WGS84' in this case.  I doubt this has anything to do 
with the problem though.  I presume it is using WGS84 internally as a default. 
 
Best regards,

Change History (2)

comment:1 by warmerdam, 19 years ago

Tyler, 

OK, a few issues coming into play.  First there was a bug in 
GDALSuggestedWarpOutput() that was causing a failure when creating the output
file for me due to the fact that corner coordinates didn't transform.  I
have committed a fix for this, though it did not seem to affect you.

Then I started getting error "geocentric transformation missing z or ellps"
instead of the "tolerance condition error" you got.  I'm not exactly sure why,
though I think I did recently make some PROJ corrections to return the 
geocentric error properly so it might be because of that, that things differ.
I determined this was in some fashion related to the fact that the Ortho
correction was setting "es" to zero, I think due to the fact that the ortho
projection only works in spherical mode. 

So I tried doing everything directly in spherical mode to avoid having 
PROJ.4 try to go through the geocentric computation to go from WGS84 to 
spherical coordinates.  This looked like:

warmerda@gdal2200[840]% gdalwarp -s_srs '+proj=latlong +ellps=sphere' -t_srs
'+proj=ortho +ellps=sphere +lat_0=60 +lon_0=-95' globe.tif ortho.tif
ERROR 1: Too many points (76 out of 84) failed to transform,
unable to compute output bounds.

Note the use of +ellps=sphere for both coordinate systems. 

I discovered that the error message is because while GDALWarpSuggestedOutput()
(the code that sets an output file size) will fallback to using an internal 
grid of points when trying to transform the region, the lower level warper
function ComputeSourceWindow() does not do this.  It turns out I have provided
specialized warp options to force this though, the SAMPLE_GRID=YES option. 
I used then and got a decent result, but with some cracking at the
anti-meridian.  This was because the sampling grid wasn't very fine so we only
realize we need data near the anti-meridian but not all the way to it.  

The solution to this was to use the SOURCE_EXTRA warp option to grab a bit
bigger window of source data for the warping.  The final result works perfectly
but it was alot of digging to figure out how to get it. 

gdalwarp -wo SOURCE_EXTRA=100 -wo SAMPLE_GRID=YES -s_srs '+proj=latlong
+ellps=sphere' -t_srs '+proj=ortho +ellps=sphere +lat_0=60 +lon_0=-95' globe.tif
ortho.tif

So, please update your GDAL from CVS or a nightly snapshot, and try with
the above options. 

I have a few things I would like to correct in the code:

 1) Figure out why the geocentric conversion failed when es was zero. 
 2) Make ComputeSourceWindow() automatically fallback to using a grid of points
    if the edge points don't all transform successfully.
 3) Automatically force a modest SOURCE_EXTRA boundary when some grid points
    don't transform to avoid cracking due to sampling error in picking the
    input window.  

Once this is all done your original transformation should work, though the
use of ellps=sphere might still be preferable to avoid any odd effects from
the internal conversion between ellipses. 


comment:2 by warmerdam, 19 years ago

I have changed gdalwarpoperation.cpp to fallback to using a grid for
sampling the source window to load if any edge pixels of the edges fail to
transform. 

I have changed gdalwarpoperation.cpp to add 10 pixels around the loaded window
if any pixels failed to transform (and SOURCE_EXTRA is not explicitly set). 

In PROJ.4 I found the problem with geocentric conversion was that 
pj_geodetic_to_geocentric was trying to transform points to geocentric even 
if they had failed the previous projection step.  Added check for HUGE_VAL 
to skip them.  The problem wasn't really related to spherical coordinates, it
is just anything that triggered a change of ellipsoid caused the geocentric
transformation to be triggered. 

All known issues now resolved in the code (GDAL and PROJ.4). 

Whew!

Note: See TracTickets for help on using tickets.