Skip to content

Commit

Permalink
Shapefile: fix wrong SRS when reading a WGS 84 SRS with a TOWGS84[0,0…
Browse files Browse the repository at this point in the history 
…,0,0,0,0,0] (fixes #3958)
  • Loading branch information
@rouault
rouault committed Jun 10, 2021
1 parent f21b6a6 commit 4eedd69
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 62 deletions.
47 changes: 12 additions & 35 deletions autotest/ogr/ogr_shape.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3808,47 +3808,24 @@ def test_ogr_shape_98():
assert ref_shx == got_shx, 'Rebuilt shx is different from original shx.'

###############################################################################
# Import TOWGS84 from EPSG when possible (#6485)
# Test WGS 84 with a TOWGS84[0,0,0,0,0,0]


def test_ogr_shape_99():
def test_ogr_shape_wgs84_with_zero_TOWGS84():

ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource('/vsimem/ogr_shape_99.shp')
lyr = ds.CreateLayer('ogr_shape_99')
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource('/vsimem/test_ogr_shape_wgs84_with_zero_TOWGS84.shp')
lyr = ds.CreateLayer('test_ogr_shape_wgs84_with_zero_TOWGS84')
ds = None
gdal.FileFromMemBuffer('/vsimem/ogr_shape_99.prj', """PROJCS["CH1903_LV03",GEOGCS["GCS_CH1903",DATUM["D_CH1903",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["False_Easting",600000.0],PARAMETER["False_Northing",200000.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Azimuth",90.0],PARAMETER["Longitude_Of_Center",7.439583333333333],PARAMETER["Latitude_Of_Center",46.95240555555556],UNIT["Meter",1.0],AUTHORITY["EPSG",21781]]""")
ds = ogr.Open('/vsimem/ogr_shape_99.shp')
gdal.FileFromMemBuffer('/vsimem/test_ogr_shape_wgs84_with_zero_TOWGS84.prj', """GEOGCS["WGS84 Coordinate System",DATUM["WGS 1984",SPHEROID["WGS 1984",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"],AUTHORITY["SBMG","LAT-LONG,LAT-LONG,WGS84,METERS"]]""")
ds = ogr.Open('/vsimem/test_ogr_shape_wgs84_with_zero_TOWGS84.shp')
lyr = ds.GetLayer(0)
got_wkt = lyr.GetSpatialRef().ExportToPrettyWkt()
expected_wkt = """PROJCS["CH1903 / LV03",
GEOGCS["CH1903",
DATUM["CH1903",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
AUTHORITY["EPSG","6149"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4149"]],
PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
PARAMETER["latitude_of_center",46.9524055555556],
PARAMETER["longitude_of_center",7.43958333333333],
PARAMETER["azimuth",90],
PARAMETER["rectified_grid_angle",90],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",600000],
PARAMETER["false_northing",200000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","21781"]]"""
ds = None

assert got_wkt == expected_wkt, ('Projections differ: got %s' % got_wkt)

ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('/vsimem/ogr_shape_99.shp')
ds = None

assert got_wkt.startswith('GEOGCS[')
assert '4326' in got_wkt

ogr.GetDriverByName('ESRI Shapefile').DeleteDataSource('/vsimem/test_ogr_shape_wgs84_with_zero_TOWGS84.shp')

###############################################################################
# Test REPACK with both implementations
Expand Down
46 changes: 19 additions & 27 deletions gdal/ogr/ogrsf_frmts/shape/ogrshapelayer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2180,33 +2180,6 @@ OGRSpatialReference *OGRShapeGeomFieldDefn::GetSpatialRef() const
}
CSLDestroy( papszLines );

// Some new? shapefiles have EPSG authority nodes (#6485) Use
// them to 'import' TOWGS84 from EPSG definition, if no
// TOWGS84 is present in the .prj (which should be the case).
// We could potentially import more, or just replace the
// entire definition
const char* pszAuthorityName = nullptr;
const char* pszAuthorityCode = nullptr;
double adfTOWGS84[7] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
if( poSRS != nullptr &&
poSRS->GetTOWGS84(adfTOWGS84, 7) == OGRERR_FAILURE &&
(pszAuthorityName = poSRS->GetAuthorityName(nullptr)) != nullptr &&
EQUAL(pszAuthorityName, "EPSG") &&
(pszAuthorityCode = poSRS->GetAuthorityCode(nullptr)) != nullptr )
{
const int nEPSGCode = atoi(pszAuthorityCode);
OGRSpatialReference oSRS;
if( oSRS.importFromEPSG(nEPSGCode) == OGRERR_NONE &&
oSRS.GetTOWGS84(adfTOWGS84, 7) == OGRERR_NONE )
{
CPLDebug(
"Shape", "Importing TOWGS84 node from EPSG definition");
poSRS->SetTOWGS84(adfTOWGS84[0], adfTOWGS84[1], adfTOWGS84[2],
adfTOWGS84[3], adfTOWGS84[4], adfTOWGS84[5],
adfTOWGS84[6]);
}
}

if( poSRS )
{
if( CPLTestBool(CPLGetConfigOption("USE_OSR_FIND_MATCHES", "YES")) )
Expand All @@ -2217,10 +2190,29 @@ OGRSpatialReference *OGRShapeGeomFieldDefn::GetSpatialRef() const
poSRS->FindMatches(nullptr, &nEntries, &panConfidence);
if( nEntries == 1 && panConfidence[0] >= 90 )
{
std::vector<double> adfTOWGS84(7);
if( poSRS->GetTOWGS84(&adfTOWGS84[0], 7) != OGRERR_NONE )
{
adfTOWGS84.clear();
}

poSRS->Release();
poSRS = reinterpret_cast<OGRSpatialReference*>(pahSRS[0]);
poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
CPLFree(pahSRS);

// If the SRS is EPSG:4326 with TOWGS84[0,0,0,0,0,0], then
// just use plain EPSG:4326
const char* pszAuthorityName = nullptr;
const char* pszAuthorityCode = nullptr;
if( adfTOWGS84 == std::vector<double>(7) &&
(pszAuthorityName = poSRS->GetAuthorityName(nullptr)) != nullptr &&
EQUAL(pszAuthorityName, "EPSG") &&
(pszAuthorityCode = poSRS->GetAuthorityCode(nullptr)) != nullptr &&
EQUAL(pszAuthorityCode, "4326") )
{
poSRS->importFromEPSG(4326);
}
}
else
{
Expand Down

0 comments on commit 4eedd69

Please sign in to comment.