[Gdal-dev] How do I make an OGRSpatialReference from this?

[Vincent Schut]

See comments inline

-------- Original Message --------
From: Christopher Barker <Chris.Barker at noaa.gov>
To: gdal-dev at lists.maptools.org
Subject: [Gdal-dev] How do I make an OGRSpatialReference from this?
Date: 07/18/2006 12:09 AM
> Hi all,
>
> In the long term, I need to be able to read a variety of
> geo-referenced raster maps, and be able to map lat/long to pixel and
> back again. I've gone through osr_tutorial.html, and been playing
> around with GDAL, all with Python. I can read the images and display
> them and I can get the sample code in the osr tutorial to work.
> However, for my two test data sets, I can't figure out how to set up
> an OGRSpatialReference.
>
> My first data set is a NOAA Raster Nautical chart (BSB). when I run
> gdalinfo on it, I get:
>
> 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)
> GCP[  2]: Id=GCP_3, Info=
>           (10587,483) -> (-88.7,29.3502,0)
> .
> .
> .
> <and a bunch more ...>
> .
> .
> .
> 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)
> Band 1 Block=10877x1 Type=Byte, ColorInterp=Palette
>   Color Table (RGB with 12 entries)
>     0: 0,0,0,255
> .
> .
> .
> It seems that the WKT is defining the datums (data?), but not the
> projection. How do i set the projection? Do I need to calculate it
> from the 50! GCPs?
This file is unprojected, but the georeferencing is defined by those
gcp's. Each gcp is a known georeferenced coordinate in a projected CRS
for a certain pixel. Using these gcp's, you can warp your file into a
projected raster. You can use gdalwarp for that, it comes with gdal. See
the help page on gdalwarp for more info on how to use it. But probably
it's just as simple as 'gdalwarp infile outfile' as gdalwarp will
recognize your input file's gcp's and their projection, and probably
(I'm not sure about this one) default to warping it to the projection
for which the gcps are defined. You can of course warp to another
projection. Depending on the type of file you have you might need to
specify a different order of warping, or maybe the special rubber sheet
warping which is especially good for gcp-georeferenced files sometimes.
Some trial and error might be needed.
The resulting output file will have a projection and correct corner
coordinates if you run gdalinfo on it.
>
> My other example is a geotiff, dumped out of ESRI ArcMap. Here's what
> I get for that:
>
> Driver: GTiff/GeoTIFF
> Size is 644, 526
> Coordinate System is `'
> Origin = (-95.395352,29.952230)
> Pixel Size = (0.00080067,-0.00080067)
> Metadata:
>   TIFFTAG_XRESOLUTION=96
>   TIFFTAG_YRESOLUTION=96
>   TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
> Corner Coordinates:
> Upper Left  ( -95.3953516,  29.9522296)
> Lower Left  ( -95.3953516,  29.5310747)
> Upper Right ( -94.8797171,  29.9522296)
> Lower Right ( -94.8797171,  29.5310747)
> Center      ( -95.1375343,  29.7416522)
> Band 1 Block=644x8 Type=Byte, ColorInterp=Red
> Band 2 Block=644x8 Type=Byte, ColorInterp=Green
> Band 3 Block=644x8 Type=Byte, ColorInterp=Blue
>
>
> For this one, I don't have a WKT, or anything else I can interpret as
> a datum or projection. I could just interpolate from the corners -- is
> that what I'm supposed to do? Is there a way to create a
> SpatialReference from this data?
So this file has some coordinates, but lacks projection info. You cannot
create projection information out of the blue, but we could do an
estimated guess. Looking at the coordinates I'd say it is in deegrees,
not in meters. If you know what area your file should contain, you can
check this. Then I'd try a simple unprojected geographic 'projection',
which is actually a simple linear mapping from lon/lat values to pixels.
I think you can add projection info using gdal_translate (the -a_srs
option). Easiest is to use 'EPSG:4326' as srs, which is a shortcut for
geographic latlong. so: 'gdal_translate -a_srs "EPSG:4326" infile
outfile'. This will give you a new file with the projection info added.
You can then check using a GIS or other reference data if the output
file's projection matches the real world situation.

Hope this helps,

Vincent.
>
> As is probably clear, I'm very new to all this!
>
> Thanks,
> -Chris
>
>
>
>
>
>