Ticket #481 (assigned defect)

Opened 4 years ago

Last modified 7 months ago

ogr2ogr: SHAPE-> MapInfo conversion - proj info get's lost

Reported by: neteler Assigned to: warmerdam (accepted)
Priority: highest Milestone:
Component: OGR_SF Version: unspecified
Severity: normal Keywords:
Cc: mloskot

Description (Last modified by mloskot)

Frank,

by chance I found a problem with projection definitions
when creating MapInfo format from SHAPE. The projection info
seems to get lost.

Here a test:

#Original file with no projection:
ogrinfo -summary ctr10000.shp ctr10000
Had to open data source read-only.
INFO: Open of `ctr10000.shp'
using driver `ESRI Shapefile' successful.

Layer name: ctr10000
Geometry: Polygon
Feature Count: 586
Extent: (1623653.000000, 4960499.000000) - (1824909.000000, 5178586.000000)
Layer SRS WKT:
(unknown)
A_CODICE: Integer (6.0)
B_NOME: String (40.0)
AGG: Integer (4.0)
DITTA: String (40.0)
DXF: String (16.0)
RVE: String (16.0)
SHP: String (16.0)
ASC: String (16.0)
VECTOR_TYP: String (16.0)

#######
# This map I wanted to fix by applying the related EPSG code. Output
#desired in MapInfo format. Using ogr2ogr to apply Monte Mario Italy 2:
ogr2ogr -a_srs '+init=epsg:26592' -f "MapInfo File" ctr10000 ctr10000.shp

#verification:
ogrinfo -summary ctr10000/ctr10000.tab ctr10000
Had to open data source read-only.
INFO: Open of `ctr10000/ctr10000.tab'
using driver `MapInfo File' successful.

Layer name: ctr10000
Geometry: Polygon
Feature Count: 586
Extent: (1623653.010000, 4960498.995000) - (1824909.000000, 5178586.005000)
Layer SRS WKT:
PROJCS["unnamed",
    GEOGCS["unnamed",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",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",2520000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1.0]]
A_CODICE: Integer (0.0)
B_NOME: String (40.0)
AGG: Integer (0.0)
DITTA: String (40.0)
DXF: String (16.0)
RVE: String (16.0)
SHP: String (16.0)
ASC: String (16.0)
VECTOR_TYP: String (16.0)

# --> It didn't pick up the EPSG definition.

#######
# Same test with SHAPE output:

ogr2ogr -a_srs '+init=epsg:26592'  ctr10000 ctr10000.shp

#verification:
ogrinfo -summary ctr10000/ctr10000.shp ctr10000
INFO: Open of `ctr10000/ctr10000.shp'
using driver `ESRI Shapefile' successful.

Layer name: ctr10000
Geometry: Polygon
Feature Count: 586
Extent: (1623653.000000, 4960499.000000) - (1824909.000000, 5178586.000000)
Layer SRS WKT:
PROJCS["Monte Mario (Rome) / Italy zone 2",
    GEOGCS["Monte Mario (Rome)",
        DATUM["Monte_Mario_Rome",
            SPHEROID["International 1924",6378388,297,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6806"]],
        PRIMEM["Rome",12.45233333333333,
            AUTHORITY["EPSG","8906"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4806"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",2520000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26592"]]
A_CODICE: Integer (6.0)
B_NOME: String (40.0)
AGG: Integer (4.0)
DITTA: String (40.0)
DXF: String (16.0)
RVE: String (16.0)
SHP: String (16.0)
ASC: String (16.0)
VECTOR_TYP: String (16.0)

# --> This looks ok.

########
#third test: convert fixed SHAPE to MapInfo:
cd ctr10000

#conversion of projected SHAPE to MapInfo:
ogr2ogr -f "MapInfo File" ctr10000 ctr10000.shp

ogrinfo -summary ctr10000/ctr10000.tab ctr10000
Had to open data source read-only.
INFO: Open of `ctr10000/ctr10000.tab'
using driver `MapInfo File' successful.

Layer name: ctr10000
Geometry: Polygon
Feature Count: 586
Extent: (1623653.010000, 4960498.995000) - (1824909.000000, 5178586.005000)
Layer SRS WKT:
PROJCS["unnamed",
    GEOGCS["unnamed",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",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",2520000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1.0]]
A_CODICE: Integer (0.0)
B_NOME: String (40.0)
AGG: Integer (0.0)
DITTA: String (40.0)
DXF: String (16.0)
RVE: String (16.0)
SHP: String (16.0)
ASC: String (16.0)
VECTOR_TYP: String (16.0)

# --> damaged projection

Probably I am missing something.

Best regards

 Markus Neteler

Change History

02/02/04 11:52:13 changed by warmerdam

Markus,

I have determined that the problem is in TABFile::SetSpatialRef() method
(gdal/ogr/ogrsf_frmts/mitab/mitab_spatialref.cpp), a part of MITAB (but
a portion I am responsible for).  

The issue is that if a coordinate system being written does not match up with
a known datum then it just falls back to using WGS84.  .tab files can actually
store any supported ellipsoid with arbitrary datum transformations it isn't
clear if the MapInfo software actually supports things not picked from their
pre-defined datum list properly.

The EPSG:26592 should look like this:
PROJCS["Monte Mario (Rome) / Italy zone 2",
    GEOGCS["Monte Mario (Rome)",
        DATUM["Monte_Mario_Rome",
            SPHEROID["International 1924",6378388,297,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6806"]],
        PRIMEM["Rome",12.45233333333333,
            AUTHORITY["EPSG","8906"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4806"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",2520000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26592"]]

Note that I am unaware of any support in MapInfo tab for non-greenwich 
meridians.  So even if I fix other aspects of the datum problem, the 
coordinate system as a whole will still not translate. 

03/22/04 20:20:44 changed by mark@keyhole.com

I'm experiencing a related problem.  I have an ESRI Shapefile of a parcel map in
Ohio.  When I convert it with ogr2ogr to a MapInfo File, some of the SRS seems
to get confused.  

Additionally, when I open up this new parcels.tab file in MapInfo Professional
6.5, and click on the projection button, it reports "Longitude / Latitude".

---------------------------------------

INFO: Open of `parcels.shp'
using driver `ESRI Shapefile' successful.
 
Layer name: parcels
Geometry: Polygon
Feature Count: 64811
Extent: (1759373.000000, 167452.265625) - (1900230.500000, 284090.218750)
Layer SRS WKT:
PROJCS["NAD_1983_StatePlane_Ohio_North_FIPS_3401_Feet",
    GEOGCS["GCS_North_American_1983",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",1968500.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",-82.5],
    PARAMETER["Standard_Parallel_1",40.43333333333333],
    PARAMETER["Standard_Parallel_2",41.7],
    PARAMETER["Latitude_Of_Origin",39.66666666666666],
    UNIT["Foot_US",0.3048006096012192]]

---------------------------------------

INFO: Open of `parcels.tab'
using driver `MapInfo File' successful.
 
Layer name: parcels
Geometry: Polygon
Feature Count: 64811
Extent: (1759373.010000, 167452.260000) - (1900230.510000, 284090.220000)
Layer SRS WKT:
PROJCS["unnamed",
    GEOGCS["unnamed",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 80",6378137,298.257222101],
            TOWGS84[0,0,0,-0,-0,-0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",40.43333333333333],
    PARAMETER["standard_parallel_2",41.7],
    PARAMETER["latitude_of_origin",39.66666666666666],
    PARAMETER["central_meridian",-82.5],
    PARAMETER["false_easting",6458320.416666667],
    PARAMETER["false_northing",0],
    UNIT["Meter",1.0]]

11/11/04 16:51:03 changed by jdunn803@gmail.com

I ran into a similar problem as what Mark reported on 3/22 and traced it to
SetSpatialRef in mitab_spatialref.cpp where the false eastings/northings
get written out.  What I found is that 
GetProjParm retrieves the item in UNITS (in Mark's case UNITS are feet, with
a dfLinearConv = 0.304xxx.  Well if you divide 1968500.0 by 0.304xxxx you get
6475328.947xxx but those aren't meters... so I think the logic in SetSpatialRef
is incorrect.  I would think that since you do set a LinearUnit in there, that
you would want to set the falses in "units" and not in meters anyway.

My suggested change (I have a diff if you want) is to remove all of the 
" / dfLinearConv" in SetSpatialRef as so:

--- mitab_spatialref.cpp	2004-11-11 15:06:40.181343500 -0500
+++ ../../../../../../temp/mitab-1.3.0/mitab/mitab_spatialref.cpp	2003-03-21
09:20:42.000000000 -0500
@@ -1035,8 +1035,8 @@
         parms[1] = poSpatialRef->GetProjParm(SRS_PP_LATITUDE_OF_CENTER,0.0);
         parms[2] = poSpatialRef->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
         parms[3] = poSpatialRef->GetProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
-        parms[4] = poSpatialRef->GetProjParm(SRS_PP_FALSE_EASTING,0.0);
-        parms[5] = poSpatialRef->GetProjParm(SRS_PP_FALSE_NORTHING,0.0);
+        parms[4] = poSpatialRef->GetProjParm(SRS_PP_FALSE_EASTING,0.0) /
dfLinearConv;
+        parms[5] = poSpatialRef->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) /
dfLinearConv;
     }

12/05/07 13:56:32 changed by gvanmaren

  • priority changed from high to highest.

Hi,

Same problem still exists?

I am trying to convert a projected shapefile to MapInfo?? format. It seems that the projection info gets lost.

When I create a MapInfo?? layer with OGR_DS_CreateLayer and OGRSpatialReferenceH defined (NZMG) -> the resulting Spatial ref for the MapInfo?? is geographic rather than the projected NZMG. It seems it falls back to using WGS84.

Is this still the same bug or is it not possible to define the NZMG projection?

Regards Gert

12/22/07 03:59:39 changed by mloskot

  • cc set to mloskot.
  • description changed.

12/29/07 11:17:29 changed by neteler

  • reporter changed from neteler@itc.it to neteler.

01/11/08 06:05:35 changed by asgerpetersen

I have just met a problem which shows how fragile it is to translate the datums by looking up the name. I had to translate a shape-file in Utm32 Etrs89 to MapInfo?, but the tab driver uses a synonym for ETRS89 which is Euref89 and furthermore it is misspelled "Euref_98". so I had to change the datum name to "Euref_98" to get it work correctly. Phew!

Isn't here an official naming convention for datums etc. in WKT? Looking up ETRS89 in "OpenGISĀ® Implementation Specification for Geographic information - Simple feature access - Part 1: Common architecture" annex b it is called "European Reference System 1989". But i have never seen this exact string used in WKT. Well, I'm confused!

Anyway, it would be nice if "ETRS_1989" was added to the code and "Euref_98" was either deleted or changed to whatever is correct.

Ogrinfo on the shape file that won't translate correctly returns:

INFO: Open of `kurver_1km_6151_537.shp'
      using driver `ESRI Shapefile' successful.

Layer name: kurver_1km_6151_537
Geometry: 3D Line String
Feature Count: 31
Extent: (537476.882888, 6151645.250881) - (537680.536825, 6152000.000000)
Layer SRS WKT:
PROJCS["ETRS_1989_UTM_Zone_32N",
    GEOGCS["GCS_ETRS_1989",
        DATUM["ETRS_1989",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",9.0],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]]
ID: Integer (8.0)
kote: Real (12.3)
objektkode: Real (11.0)

Ogrinfo on the tab file (which works correctly in MapInfo?) shows:

Had to open data source read-only.
INFO: Open of `kurver_1km_6151_537.tab'
      using driver `MapInfo File' successful.

Layer name: kurver_1km_6151_537
Geometry: Line String
Feature Count: 31
Extent: (537476.884067, 6151645.250403) - (537680.539929, 6151999.998195)
Layer SRS WKT:
PROJCS["unnamed",
    GEOGCS["unnamed",
        DATUM["Euref_98",
            SPHEROID["GRS 80",6378137,298.257222101],
            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",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1.0]]
ID: Integer (0.0)
kote: Real (0.0)
objektkode: Real (0.0)

Best regards Asger