[Xastir] Ordnance Survey (OSGB) "OS OpenMap Local"

[Tom Russo]

On Wed, May 02, 2018 at 09:39:37AM -0500, we recorded a bogon-computron collision of the <kg4wsv at gmail.com> flavor, containing:
> On Wed, May 2, 2018 at 9:20 AM, Tom Russo <russo at bogodyn.org> wrote:
> 
> > It would be more work than anyone's been willing to put in in the nearly
> > 14 years I've been using Xastir.
> 
> OTOH, this is how much effort it would be to convert a directory full
> of maps to the commonly accepted standard WGS84:
> 
> $ cd directory/to/maps
> $ mkdir ../new
> $ for x in *.tif; do gdalwarp -t_srs EPSG:4326 $x ../new/$x; done
> 
> which is one reason such a change request wouldn't get much traction.

What he said.

> I'm sure no one would protest if someone wanted to volunteer to do it.
> :)

Certainly not.  If a user were to pick apart the map_tif.c code and figure
out how to make it really do multiple projections/coordinate systems right,
the effort would certainly not be rejected.  But basically, this issue has
been well known for at least 14 years, and nobody's taken up the flag to 
address it directly.  Instead, we just recommend the cop-out approach of
fixing your map data for Xastir, instead of extending Xastir's compatability
with every one of the nearly infinite set of map projections and coordinate
systems.

And if you're really interested in why this is hard... get ready for a novel.

The GeoTIFF code is something that has a lot of assumptions in it based on
the original use case:  Curt wanted to be able to read USGS 7.5' quad files
for his hunting and Search and Rescue work.  There are many assumptions and
approximations in this code that make this work quite well (and 
even efficiently) for these large scale maps in UTM coordinates, but that fail 
badly for a lot of other map projections and scales.

You'll notice that in configure.ac and in much of the Xastir documentation, 
GeoTIFF support is repeatedly referred to as if it is really a US-centric 
feature for USGS data products, even though GeoTIFF is a perfectly general 
format for worldwide use.  That's why -- it was really only written for this
use case, even though sometimes it works for others.

It's not really that Xastir supports multiple coordinate systems --- the 
GeoTIFF code really doesn't know anything about projected coordinate systems,
and just uses the GeoTIFF library's function that converts the PCS coordinates 
of the corners of the image to lat/lon, and then converts those lat/long
corner points to the "xastir coordinate system" (which is really just a 
scaled representation of the lat/lon so they can be integers).  This function
is always returning good values for lat/lon for the corners.  It's just that
Xastir's assumptions don't always work after that.

After obtaining the lat/lon values of the corners from the GeoTIFF library,
Xastir then uses these corner coordinates and the image size in pixels to 
work out where to put each pixel on the screen.

The core assumption that tends to break things is that a scan line of the image 
corresponds to a straight line (rotated and stretched, but straight) in  the 
Xastir coordinate system.  Each pixel in the scan line is simply plopped down
along this line as a rectangle of appropriate size (ideally one pixel, but
when zoomed in close, an actual rectangle).

This assumption is fine for large scale maps like 7.5' quads in UTM --- at
this scale, a scan line of the image is very, very close to a straight line
in lat/long --- but very wrong for small scale maps (i.e. those that cover 
large regions) or for other projections, where a line in the image is actually 
a curve in lat/lon.  The result is usually that the map corners are placed 
somewhat correctly and the innards of the map are misplaced.  This was 
blindingly obvious the first time someone tried to read in GeoTIFF files for 
FAA sectional aviation charts which are very small scale compared to a USGS 
7.5' quad --- they'd display, but they were wrong.

But it could also be that the assumption is so bad for how the Ordnance Survey
maps are made that it can't even get the corners right, leading to the sort of
error you have here.  

 Aside added late in writing/editing this email:

  Actually, there's a second assumption I've just spotted by reading the 
  map_tif code that is also biting here.  There is an assumption (grossly 
  incorrect in this case) that if the image is using any datum other than 
  WGS84 or NAD83 that it must in fact be in NAD27 and must have a datum shift 
  applied to it.  This is very much a USGS 7.5' quadrangle assumption, and is 
  really invalid here.  Your map's datum is OSGB 1936.  

  I have no idea what the home-grown datum shifting code in Xastir is doing to 
  your coordinates when it assumes they are NAD27 and applies an approximate 
  datum shift to them.  NAD27 is only valid in North America, and the 
  approximate Molodensky transformation that is used by Xastir (instead of the
  more correct NADCON datum shifting grids) will probably be *very* wrong
  in Great Britain.  This may be the *real* reason that the code is puking
  hard on the coordinates and rejecting the raster.

Because the GeoTIFF code makes this really simplifying (but bad) assumption
that it can just linearly place scan lines of the image, the easiest way to 
get maps into Xastir that are in other projections is simply to warp them into 
geographic coordinates (lat/lon) using gdalwarp --- in which case the 
assumption that a line of the image is a straight line in lat/lon is trivially 
true.   

The initial map_tif.c code couldn't work with such maps (it flagged any geotiff
without a projected coordinate system as defective), but in 2004 I hacked 
it (almost trivially) to make them work.  It was mostly just a matter of 
removing a test that refused to read a map without a projected coordinate 
system (PCS), and skipping a bit of code that transformed coordinates from the 
PCS to lat/lon (because they were already *IN* lat/lon).  There was also
a little work to be done to adapt the rendering optimization trick that was
coded for projected maps, where it skips scan lines of the map if you're zoomed
out too far for them to be distinct in the rendered image.

Making the code in map_tif.c work properly for rasters in anything other than 
the large-scale UTM (for which it was originally written), or for geographic 
coordinates (lat/lon) (for which the work is almost trivial) would be difficult.

The best path forward would actually be to use GDAL itself to read and warp 
the rasters instead of the hand-coded and optimized approximations that are in
map_tif.c.  Unfortunately, there was a tiny amount of work back in 2003 to 
integrate GDAL and OGR into Xastir, but that work stalled almost immediately 
for the raster side (GDAL) and only got part way for the vector side (OGR) 
before development stopped sometime in 2004.  This is why we generally tell
users not to bother building Xastir with GDAL support unless they need gdal
for something else (like, say, the various programs it provides like gdalwarp).

If someone *really* wanted to add value to Xastir, it would probably be to
do a proper implementation of raster rendering, using GDAL to read the
images in instead of picking it apart and using homegrown approximations to
render it.  This is not a small thing.  One of the reasons it never got worked 
on much back in the 2003-2004 days was that warping a single raster from one 
coordinate system to another is really processor-intensive work, and Xastir 
raster rendering on 2003-2004 vintage computers was already insanely slow.  

When you get around to running the gdalwarp commands suggested in this thread,
even on today's computers it's not exactly blindingly fast.  Proper warping
of raster images (especially high resolution images) from one coordinate system
to another is still time consuming.  It is better just to do it as a 
preprocessing step and be done with it.

-- 
Tom Russo    KM5VY
Tijeras, NM  

 echo "prpv_a'rfg_cnf_har_cvcr" | sed -e 's/_/ /g' | tr [a-m][n-z] [n-z][a-m]