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
_______________________________________________
grassuser mailing list
[email protected]
http://grass.itc.it/mailman/listinfo/grassuser
