RMF: Incorrect zone number detection

[tsarioann]

Currently in frmts/rmf/rmfdataset.cpp zone number is detected from x-coordinate of left border. On practice left border sometimes is in another zone. This patch is for detecting zone number from center-coordinate:

--- frmts/rmf/rmfdataset.cpp	2012-05-17 07:30:55.000000000 +0400
+++ frmts/rmf/rmfdataset.cpp	2012-07-20 18:44:07.264162201 +0400
@@ -1486,8 +1486,8 @@
         {
             if ( poDS->sExtHeader.nZone == 0 )
             {
-                padfPrjParams[7] =
-                    floor((poDS->sHeader.dfLLX - 500000.0 ) / 1000000.0);
+                double centerXCoord = poDS->sHeader.dfLLX + (poDS->nRasterXSize * poDS->sHeader.dfPixelSize / 2.0);
+                padfPrjParams[7] = floor((centerXCoord - 500000.0 ) / 1000000.0);
             }
             else
                 padfPrjParams[7] = poDS->sExtHeader.nZone;

P.S. Can't commit myself, cause I don't have permissions.

[Even Rouault]

The idea behind the patch and the implementation seem reasonable, but not being a specialist of the RMF format/driver, I'm not sure if it's "universally" correct. CC'ing Andrey Kiselev who is the author of the driver if he has some though about that ? Unless that you (tsarioann) have good knowledge of the format ?

[tsarioann]

Replying to rouault:

The idea behind the patch and the implementation seem reasonable, but not being a specialist of the RMF format/driver, I'm not sure if it's "universally" correct. CC'ing Andrey Kiselev who is the author of the driver if he has some though about that ? Unless that you (tsarioann) have good knowledge of the format ?

I don't have good knowledge of the format. But Panorama loads my rsw files (where left border is in other zone) correctly and GDAL does not. This means that Panorama and GDAL have different implementation of zone detection. I hope Andrey Kiselev has sources of Panorama and can reimplement zone detection from there to GDAL.

[Even Rouault]

I guess we should align on Panorama behaviour as being the reference.

But reviewing the patch now (and even the version before your patch), there's something strange : how can we guess a UTM zone number from the x coordinate of a point of the image ? If a point has x = 500 000, this doesn't give you much information : it might be located in any of the 60 UTM zones. Unless the x coordinates in .rsw go outside the typical range (which is something like [200 000,800 000] usually) and are counted as meters from the -180 meridian. But if that was the case, they should be corrected once the zone number is guessed, and I don't see where it is done in the code. So I'm obviously missing something here, but I fail to see what.

[tsarioann]

My rsw files use Pulkovo_1942 datum. So here are corner coordinates from two files:

'''M-38-XXII.rsw:'''
Corner Coordinates:
Upper Left  ( 8495000.000, 5546420.000) ( 44d55'48.69"E, 50d 2'57.14"N)
Lower Left  ( 8495000.000, 5461790.000) ( 44d55'52.58"E, 49d17'17.93"N)
Upper Right ( 8576700.000, 5546420.000) ( 46d 4'14.79"E, 50d 2'39.43"N)
Lower Right ( 8576700.000, 5461790.000) ( 46d 3'15.14"E, 49d17' 0.69"N)

'''M-38-XXIII.rsw'''
Corner Coordinates:
Upper Left  ( 8567680.000, 5547865.000) ( 45d56'42.45"E, 50d 3'30.13"N)
Lower Left  ( 8567680.000, 5462275.000) ( 45d55'49.19"E, 49d17'20.22"N)
Upper Right ( 8648390.000, 5547865.000) ( 47d 4'17.91"E, 50d 2'37.42"N)
Lower Right ( 8648390.000, 5462275.000) ( 47d 2'21.24"E, 49d16'28.91"N)

[Even Rouault]

And which SRS does gdalinfo report before and after your correction ?

[tsarioann]

Before:

Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Pulkovo 1942",
        DATUM["Pulkovo_1942",
            SPHEROID["Krassowsky 1940",6378245,298.3,
                AUTHORITY["EPSG","7024"]],
            TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
            AUTHORITY["EPSG","6284"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4284"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",7500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (8495000.000000000000000,5546420.000000000000000)

After:

Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Pulkovo 1942",
        DATUM["Pulkovo_1942",
            SPHEROID["Krassowsky 1940",6378245,298.3,
                AUTHORITY["EPSG","7024"]],
            TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
            AUTHORITY["EPSG","6284"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4284"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",45],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",8500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (8495000.000000000000000,5546420.000000000000000)

central_meridian and false_easting are incorrect in the first variant.

[Even Rouault]

ok, I have applied your patch in trunk (r24731) and branches/1.9 (r24732). Thanks

[tsarioann]

Thank you.