[GRASS-user] OT: transform mapinfo tab file with coordsys to world file

[Moritz Lennert]

On 27/10/07 04:45, Hamish wrote:
> Moritz Lennert wrote:
>> We have a series of georeferenced .tif files accompanied by .tab 
>> files which look like this:
>> 
>> Definition Table File "s5_15 nw3_b-iii.tif" Type "RASTER" 
>> (513881,9523020) (1008,1098) Label "Pt 1", (519881,9523020) 
>> (10493,1111) Label "Pt 2", (519881,9517020) (10451,10523) Label "Pt
>>  3", (513881,9517020) (990,10509) Label "Pt 4" CoordSys Earth 
>> Projection 8, 104, "m", 15, 0, 0.9996, 500000, 10000000 Units "m"
>> 
>> However, gdal does not support the Mapinfo's CoordSys: "only 
>> control points used, Coordsys ignored" 
>> (http://www.gdal.org/frmt_gtiff.html).
>> 
>> I guess that if we could transform these files into ESRI world 
>> files (http://www.gdal.org/frmt_various.html#WLD), we could import 
>> them easily with GDAL. But how to do this, as ESRI world files give
>> transformation parameters, not points and projection info ?
>> 
>> Or is this the wrong approach ?
>> 
>> We can obviously easily import each file into an XY location and 
>> then run i.rectify on the basis of the available points, but I was 
>> wondering whether there is a better way.
> 
> 
> 1) check if the georeferecing info is registered in the tiff file (if
>  it is a GeoTIFF) with 'gdalinfo'. If it is, just use r.in.gdal.

This is the relevant output:

Coordinate System is `'
GCP Projection =
PROJCS["unnamed",GEOGCS["unnamed",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563],
TOWGS84[0,0,0,-0,-0,-0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],
PARAMETER["false_northing",10000000],UNIT["Meter",1]]


So it actually seems to recognizes everything. But r.in.gdal does not
seem to use this info during import.

> 
> 2) If gdalinfo shows no georeferecing info in the file, you can add 
> those four ground control points with 'gdal_translate' -gcp'. Note 
> the "line" parameter (y pixel) may start from the top, not the 
> bottom, so you might need to do height-line_number. (only try that if
>  it fails)

The gcp's are imported when we import the file into an XY-location.

> 
> A quick search of the EPSG file shows 6 possibles:
[...]
> 3) Do you know what projection it is supposed to be in?


We know what the projection is. It is all contained in the MapInfo string:
Projection 8, 104, "m", 15, 0, 0.9996, 500000, 10000000 Units "m"

which translates to

PROJCS["Transverse Mercator",
     GEOGCS["wgs84",
         DATUM["WGS_1984",
             SPHEROID["wgs84",6378137,298.257223563],
             TOWGS84[0.000,0.000,0.000]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Transverse_Mercator"],
     PARAMETER["latitude_of_origin",0],
     PARAMETER["central_meridian",15],
     PARAMETER["scale_factor",0.9996],
     PARAMETER["false_easting",10000000],
     PARAMETER["false_northing",500000],
     UNIT["meters",1]]


But as can be seen above, gdal does not recognize this.

> If you can find the EPSG code for that you can add that into the
> GeoTIFF with 'gdal_translate -a_srs'

Using -a_srs with the above WKT definition in a file (gdal_translate 
-a_srs proj S5_15\ NW3_B-IV.tif brol.tif), I now get:

Coordinate System is `'
GCP Projection = PROJCS["Transverse Mercator",GEOGCS["WGS 
84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.2572235629972,AUTHORITY["EPSG","7030"]], 
AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433], 
AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0], 
PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",10000000], 
PARAMETER["false_northing",500000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

when running gdalinfo. So, the GCP projection info is more complete, but 
the file is still not recognized as georeferenced when importing with 
'r.in.gdal location='

Maybe I'm doing something wrong...

> or create a new GRASS location using those parameters. Then import 
> with 'r.in.gdal -o'

This does not work, the map is imported in its xy-extents, not the
projected extents.

> 
> 
> Seeing you do have some numbers there already, importing to a XY 
> location and using i.rectify should be a last resort. (note, you can 
> use those four given tie points in the i.rectify POINTS file, it's 
> plain text edit it by hand after setting a few rough points to see 
> the format/order)

As mentioned above, the existing GCP's are imported by r.in.gdal.

So, we can automate r.in.gdal & i.rectify, but I was wondering whether 
there was a better way...

Moritz