[gdal-dev] Puzzling results from gdalinfo with BSB nautical charts

[Even Rouault]

Jack,

your understanding of the BSB driver source code is correct. It will always 
return WGS84 as the GCP projection.
NAD83 and WGS84 are generally considered as equivalent/compatible, so this 
should be OK for your use case.

But I've no idea if there are BSB datasets that are not NAD83 or WGS84. In 
which case, the BSB driver would have to be improved to parse the appropriate 
line.

Best regards,
Even

Le Tuesday 30 December 2008 05:53:10 Jack Kirkpatrick, vous avez écrit :
> Greetings,
>
> I have been exploring the GDAL utilities with a view to creating a general
> BSB to RLE bitmap converter that also generates a separate WGS84 worldfile,
> appropriately transformed for the datum etc of the original BSB file.
>
> I believe I am getting some puzzling results from the gdalinfo utility that
> I will try to explain as follows:
>
> --------------------------------------
>
> (1) Environment is Windows XP SP3 and I have installed FWTools 2.2.9.
>
>
> --------------------------------------
>
> (2) I have obtained NOAA BSB chart no 11006 (KEY WEST TO THE MISSISSIPPI
> RIVER) from:
>
> http://www.nauticalcharts.noaa.gov/mcd/Raster/download_agreement.htm
>
> scroll to the bottom then click [Proceed to NOAA RNC(R) (BSB) Downloads]
> then click [Textual], the rest is pretty obvious.
>
>
> --------------------------------------
>
> (3) After downloading and unzipping to [C:\NOAA BSB\11006], I inspect the
> file [C:\NOAA BSB\11006\11006_1.KAP] with the Windows Wordpad applet, the
> first several lines read as follows:
>
> !
> CRR/
> VER/3.0
> BSB/NA=KEY WEST TO THE MISSISSIPPI RIVER,NU=44,RA=10967,9389,DU=254
> KNP/SC=875000,GD=NAD83,PR=MERCATOR,PP=26.127,PI=60.000,SP=,SK=0.0000000
> TA=90.0000000,UN=FATHOMS,SD=MEAN LOWER LOW WATER,DX=87.50,DY=87.50
> ...
>
> I'm guessing "GD=NAD83,PR=MERCATOR" means the file contents are to NAD83
> datum in a Mercator projection.
>
> This file also contains 92 ground control points of the form:
>
> REF/1,275,8522,24.1668694444,-90.0831111111
> REF/2,275,300,30.5834527778,-90.0831111111
> REF/3,2427,300,30.5834527778,-88.1999805556
> ...
> REF/92,275,1069,30.0000583333,-90.0831111111
>
> --------------------------------------
>
> (4) I launch the FWTools shell via the [FWTools Shell] shortcut on the
> desktop and within the FWTools shell I key:
>
> gdalinfo -nomd -noct "C:\NOAA BSB\11006\11006_1.kap" > "C:\NOAA
> BSB\11006\11006_1.txt"
>
>
> --------------------------------------
>
> (5) I inspect the resulting file [C:\NOAA BSB\11006\11006_1.txt] which
> shows:
>
> Driver: BSB/Maptech BSB Nautical Charts
> Files: C:\NOAA BSB\11006\11006_1.kap
> Size is 10967, 9389
> Coordinate System is `'
> GCP Projection = GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.257223563,AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AU
>THORITY["EPSG",6326]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["deg
>ree",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long
>",EAST],AUTHORITY["EPSG",4326]] GCP[ 0]: Id=GCP_1, Info=
> (275,8522) -> (-90.0831111111,24.1668694444,0)
> GCP[ 1]: Id=GCP_2, Info=
> (275,300) -> (-90.0831111111,30.5834527778,0)
> GCP[ 2]: Id=GCP_3, Info=
> (2427,300) -> (-88.1999805556,30.5834527778,0)
> ...
> GCP[ 91]: Id=GCP_92, Info=
> (275,1069) -> (-90.0831111111,30.0000583333,0)
> Corner Coordinates:
> Upper Left ( 0.0, 0.0)
> Lower Left ( 0.0, 9389.0)
> Upper Right (10967.0, 0.0)
> Lower Right (10967.0, 9389.0)
> Center ( 5483.5, 4694.5)
> Band 1 Block=10967x1 Type=Byte, ColorInterp=Palette
> Color Table (RGB with 12 entries)
>
> <<<<< i.e. somehow gdalinfo has decided that the files contents are to
> WGS84 datum >>>>>
>
> --------------------------------------
>
> (6) I rummaged through the gdal source and found
> [C:\gdal-1.5.3\frmts\bsb\...] which appears to contain the source files for
> the GDAL BSB driver(?).
>
> Of particular interest is the file [C:\gdal-1.5.3\frmts\bsb\bsbdataset.cpp]
> which contains the line:
>
> pszGCPProjection = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS
> 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",7030]],TOWGS84[0,0,0,0,0,0,0]
>,AUTHORITY[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],
>UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],AXIS[\"Lat\",NO
>RTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",4326]]";
>
> which is a dead ringer for the [GSP Projection =...] produced by gdalinfo.
>
> I am not a C programmer but reading
> [C:\gdal-1.5.3\frmts\bsb\bsbdataset.cpp] and
> [C:\gdal-1.5.3\frmts\bsb\bsb_read.c] it looks to me suspiciously like the
> GDAL BSB driver forces the gdalinfo datum to WGS84 regardless of what it
> truly is.
>
> --------------------------------------
>
> (7) I am led to believe that NAD83 is to all intents and purposes WGS84 so
> I guess there is no harm done if all BSB charts were NAD83 or similarly
> "WGS84 compatible".
>
> However if you go to:
>
> http://www.digitalocean.ca/digitaloceane-store/bsbraster.asp
>
> you can obtain BSB charts for areas totally divorced from NOAA (and
> therefor presumably NAD83 and similar).
>
> My concern is that I could be in the situation of trying to use the GDAL
> utilities to convert (say) some New Zealand BSB charts that use a datum
> that is not "WGS84 compatible" but the GDAL utilities force the datum to
> WGS84 without any appropriate transformation.
>
> --------------------------------------
>
> Any and all comments welcomed - I do hope I haven't missed something
> embarrassingly simple in this great monologue...