[Gdal-dev] How do i create a Mercator Projection from a BSB file?

[Christopher Barker]

Hi all,

I'm trying to figure out how to map from pixel coords to lat-long (and 
back again) from a noaa raster nautical chart (BSB). gdalinfo gives this 
data that's relevant:

Driver: BSB/Maptech BSB Nautical Charts
Size is 10877, 9391
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],AUTHORITY["EPSG",6326]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",4326]]
GCP[  0]: Id=GCP_1, Info=
           (248,8558) -> (-89.55,28.7669,0)
GCP[  1]: Id=GCP_2, Info=
           (248,483) -> (-89.55,29.3502,0)
< and a bunch more GCPs >

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 9391.0)
Upper Right (10877.0,    0.0)
Lower Right (10877.0, 9391.0)
Center      ( 5438.5, 4695.5)


As far as I can tell, there is only Geographic Coordinate System info in 
the WKT. Here it is expanded for easier reading:
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",NORTH],
     AXIS["Long",EAST],
     AUTHORITY["EPSG",4326]]

There doesn't seem to be any projection info anywhere. Except maybe 
where it says:

UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",9108]]

What does the number after "degree" mean?

However, there is more data in the top of the *.kap file, that the 
driver doesn't seem to extract:

PTC/Maptech,KAPLib V5.0.0.4,12/20/2005
BSB/NA=MISSISSIPPI RIVER DELTA,NU=57,RA=10877,9391,DU=254
KNP/SC=80000,GD=NAD83,PR=MERCATOR,PP=29.133,PI=10.000,SP=Louisiana 
South,40000
     SK=0.0000000,TA=90.0000000,UN=FEET,SD=MEAN LOWER LOW 
WATER,DX=8.00,DY=8.00
KNQ/EC=RF,GD=NARC,VC=UNKNOWN,SC=MLLW,PC=MC,P1=UNKNOWN,P2=29.133
     P3=NOT_APPLICABLE,P4=NOT_APPLICABLE,GC=KK,RM=POLYNOMIAL
CED/SE=72,RE=1,ED=10/01/2004
NTM/NE=72.63,ND=12/17/2005,BF=OFF,BD=10/01/2004

Here's what I recognise and think is relevant from that:

Scale:80000
Projection: PR=MERCATOR
Reference Latitude: 29.133

Now I've dug into the osr.py, and found that SetMercator requires:

SetMercator( self, clat, clong, scale, fe, fn )

Which I think means:

#  which I think means;
#  	clat --- Latitude of the projection origin.
#       clong --- Longitude of the projection origin (central meridian?)
#	scale --- Scale factor applies to the projection origin.
#	fe --- False easting.
#	fn --- False northing."""

So where to I get that data from the above?

I think I got clat from there: 29.133
clong, I have no idea.

Scale is given as 80000, and that's the scale of the paper cahrt, but I 
don't know how to convert that to pixel coords.

Is there a False easting or False northing for this chart?

Any help would be appreciated.

Also, it seems that the BSB driver could be extracting more info than it 
does -- is there a reason for that? or does someone just need to write 
the code? If the later, than if someone can tell me what fields in 
GDALDataset should be filled with what from the *.kap[ file, I could 
probably submit a patch.

By the way, is there an easier way to do this? It seems a common need to 
create an OGRSpatialReference from the data in a GDALDataset. Is there 
an existing method to do that?

thanks,

-Chris


-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov