289 строки
11 KiB
Python
Исполняемый файл
289 строки
11 KiB
Python
Исполняемый файл
#!/usr/bin/env pytest
|
|
###############################################################################
|
|
# $Id$
|
|
#
|
|
# Project: GDAL/OGR Test Suite
|
|
# Purpose: Test read/write round-tripping of SRS for HFA/Imagine format.
|
|
# Author: Even Rouault, <even dot rouault at spatialys.com>
|
|
#
|
|
###############################################################################
|
|
# Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
|
|
#
|
|
# Permission is hereby granted, free of charge, to any person obtaining a
|
|
# copy of this software and associated documentation files (the "Software"),
|
|
# to deal in the Software without restriction, including without limitation
|
|
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
|
# and/or sell copies of the Software, and to permit persons to whom the
|
|
# Software is furnished to do so, subject to the following conditions:
|
|
#
|
|
# The above copyright notice and this permission notice shall be included
|
|
# in all copies or substantial portions of the Software.
|
|
#
|
|
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
|
|
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
# DEALINGS IN THE SOFTWARE.
|
|
###############################################################################
|
|
|
|
import gdaltest
|
|
import pytest
|
|
|
|
from osgeo import gdal, osr
|
|
|
|
###############################################################################
|
|
# Write a HFA/Imagine and read it back to check its SRS
|
|
crs_list = [
|
|
[3814, False, False], # tmerc
|
|
[
|
|
28991,
|
|
False,
|
|
True,
|
|
], # sterea / Oblique Stereographic # requires ESRI PE as we have no mapping for it
|
|
# [2046, False, False], # tmerc south oriented DISABLED. Not sure about the axis
|
|
[3031, False, False], # stere, "Polar Stereographic (variant B)"
|
|
[
|
|
5041,
|
|
True,
|
|
False,
|
|
], # stere, "Polar Stereographic (variant A)". failure_expected because we ignore the scale factor
|
|
[6931, False, False], # laea
|
|
[
|
|
2062,
|
|
False,
|
|
True,
|
|
], # lcc "Lambert Conic Conformal (1SP)". We morph it to LCC_2SP in the Imagine representation
|
|
[3943, False, False], # lcc "Lambert Conic Conformal (2SP)"
|
|
# [2065, True, False], # krovak South-West
|
|
[5221, True, False], # krovak east-north
|
|
[2066, False, True], # cass # issue with Clarke's link unit
|
|
[2964, False, False], # aea
|
|
[3410, False, False], # cea
|
|
[3786, False, False], # eqc
|
|
[3000, False, False], # merc, "Mercator (variant A)" with scale_factor!=1
|
|
[3832, False, False], # merc, "Mercator (variant A)" with scale_factor==1
|
|
[
|
|
5641,
|
|
False,
|
|
True,
|
|
], # merc, "Mercator (variant B)". We morph it to Mercator Variant A in the Imagine representation
|
|
[27200, False, False], # nzmg
|
|
[6842, False, False], # omerc, "Hotine Oblique Mercator (variant A)"
|
|
[2057, False, False], # omerc, "Hotine Oblique Mercator (variant B)"
|
|
[29100, False, False], # poly
|
|
[2056, False, False], # somerc
|
|
[2027, False, False], # utm
|
|
[4326, False, False], # longlat
|
|
["ESRI:54049", False, False], # Vertical Perspective
|
|
[5472, False, False], # Polyconic
|
|
["ESRI:102010", False, False], # eqdc
|
|
["ESRI:102237", False, False], # aeqd
|
|
["ESRI:102034", False, False], # Gnomonic
|
|
[
|
|
"ESRI:102035",
|
|
True,
|
|
False,
|
|
], # Orthographic # failure_expected because we return Orthographic whereas input is Orthographic (Spherical)
|
|
[
|
|
"+proj=ortho +lat_0=90 +lon_0=1 +x_0=2 +y_0=3 +datum=WGS84 +units=m +no_defs",
|
|
False,
|
|
False,
|
|
], # Orthographic
|
|
["ESRI:102011", False, False], # Sinusoidal
|
|
["ESRI:53003", False, False], # Miller
|
|
["ESRI:53029", False, False], # Van Der Grinten
|
|
["ESRI:53030", False, False], # Robinson
|
|
["ESRI:53009", False, False], # Mollweide
|
|
["ESRI:53010", False, False], # Sphere_Eckert_VI
|
|
["ESRI:53011", False, False], # Sphere_Eckert_V
|
|
["ESRI:53012", False, False], # Sphere_Eckert_IV
|
|
["ESRI:53013", False, False], # Sphere_Eckert_III
|
|
["ESRI:53014", False, False], # Sphere_Eckert_II
|
|
["ESRI:53015", False, False], # Sphere_Eckert_I
|
|
["ESRI:53016", False, False], # Sphere_Gall_Stereographic
|
|
[28191, False, False], # Cassini Soldner
|
|
["ESRI:53031", False, False], # Two_Point_Equidistant
|
|
["ESRI:102163", False, False], # Bonne
|
|
["ESRI:53023", False, False], # Sphere_Loximuthal
|
|
["ESRI:53022", False, False], # Sphere_Quartic_Authalic
|
|
["ESRI:53018", False, False], # Sphere_Winkel_I
|
|
["ESRI:53019", False, False], # Sphere_Winkel_II
|
|
[
|
|
"ESRI:102421",
|
|
True,
|
|
False,
|
|
], # Equidistant Cylindrical failure_expected because we return Equidistant Cylindrical whereas input is Equidistant Cylindrical (Spherical)
|
|
[
|
|
"+proj=eqc +lat_ts=22.94791772 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs",
|
|
False,
|
|
False,
|
|
], # Same as above but really with Equidistant Cylindrical
|
|
["ESRI:53043", False, False], # Sphere_Aitoff
|
|
["ESRI:53046", False, False], # Sphere_Craster_Parabolic
|
|
["ESRI:53045", False, False], # Flat_Polar_Quartic
|
|
["ESRI:53048", False, False], # Sphere_Times
|
|
["ESRI:53042", False, False], # Sphere_Winkel_Tripel_NGS
|
|
["ESRI:53044", False, False], # Sphere_Hammer_Aitoff
|
|
[
|
|
"ESRI:53025",
|
|
False,
|
|
False,
|
|
], # Sphere_Hotine, Hotine_Oblique_Mercator_Two_Point_Natural_Origin
|
|
]
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
"srs_def,failure_expected,read_with_pe_string",
|
|
crs_list,
|
|
ids=[str(r[0]) for r in crs_list],
|
|
)
|
|
def test_hfa_srs(srs_def, failure_expected, read_with_pe_string):
|
|
sr = osr.SpatialReference()
|
|
sr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
|
|
if isinstance(srs_def, str):
|
|
sr.SetFromUserInput(srs_def)
|
|
else:
|
|
sr.ImportFromEPSG(srs_def)
|
|
|
|
ds = gdal.GetDriverByName("HFA").Create("/vsimem/TestHFASRS.img", 1, 1)
|
|
ds.SetSpatialRef(sr)
|
|
ds = None
|
|
|
|
with gdaltest.config_options(
|
|
{} if read_with_pe_string else {"HFA_USE_ESRI_PE_STRING": "NO"}
|
|
):
|
|
ds = gdal.Open("/vsimem/TestHFASRS.img")
|
|
sr2 = ds.GetSpatialRef()
|
|
ds = None
|
|
|
|
gdal.Unlink("/vsimem/TestHFASRS.img")
|
|
|
|
# For EPSG:2065. Those 2 datums are translated into D_S_JTSK in ESRI WKT... So for the purpose of
|
|
# comparison, substitute one for another
|
|
if (
|
|
sr.ExportToWkt().find(
|
|
'"System_Jednotne_Trigonometricke_Site_Katastralni_Ferro"'
|
|
)
|
|
!= -1
|
|
and sr2.ExportToWkt().find('"System_Jednotne_Trigonometricke_Site_Katastralni"')
|
|
!= -1
|
|
):
|
|
wkt2 = sr2.ExportToWkt().replace(
|
|
'"System_Jednotne_Trigonometricke_Site_Katastralni"',
|
|
'"System_Jednotne_Trigonometricke_Site_Katastralni_Ferro"',
|
|
)
|
|
sr2.SetFromUserInput(wkt2)
|
|
|
|
if (
|
|
not isinstance(srs_def, str)
|
|
and srs_def == 4326
|
|
and sr2.GetAuthorityCode(None) != "4326"
|
|
) or sr.IsSame(sr2) != 1:
|
|
|
|
if osr.GetPROJVersionMajor() * 100 + osr.GetPROJVersionMinor() < 630:
|
|
if srs_def == "ESRI:54049":
|
|
if (
|
|
sr2.ExportToProj4()
|
|
!= "+proj=nsper +lat_0=0 +lon_0=0 +h=35800000 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
|
|
):
|
|
print(sr)
|
|
print(sr2)
|
|
print(sr2.ExportToProj4())
|
|
assert False, "did not get expected SRS"
|
|
else:
|
|
return
|
|
else:
|
|
sr.ImportFromWkt(sr.ExportToWkt())
|
|
sr2.ImportFromWkt(sr2.ExportToWkt())
|
|
|
|
if sr.IsSame(sr2) != 1:
|
|
if failure_expected:
|
|
pytest.xfail(
|
|
"did not get expected SRS. known to be broken currently. FIXME!"
|
|
)
|
|
|
|
print(sr)
|
|
print(sr2)
|
|
assert False, "did not get expected SRS"
|
|
|
|
|
|
def test_hfa_srs_wisconsin_tmerc():
|
|
|
|
ds = gdal.Open(
|
|
"data/esri_103300_NAD_1983_HARN_WISCRS_Adams_County_Meters_transverse_mercator.img"
|
|
)
|
|
wkt = ds.GetProjectionRef()
|
|
sr = osr.SpatialReference()
|
|
sr.SetFromUserInput(wkt)
|
|
assert sr.GetAuthorityCode(None) == "103300"
|
|
|
|
|
|
def test_hfa_srs_NAD83_UTM():
|
|
sr = osr.SpatialReference()
|
|
sr.ImportFromEPSG(26915)
|
|
|
|
ds = gdal.GetDriverByName("HFA").Create("/vsimem/TestHFASRS.img", 1, 1)
|
|
ds.SetProjection(sr.ExportToWkt())
|
|
ds = None
|
|
|
|
ds = gdal.Open("/vsimem/TestHFASRS.img")
|
|
wkt = ds.GetProjectionRef()
|
|
assert ds.GetSpatialRef().GetAuthorityCode(None) == "26915"
|
|
ds = None
|
|
|
|
gdal.Unlink("/vsimem/TestHFASRS.img")
|
|
|
|
assert "TOWGS84" not in wkt
|
|
|
|
|
|
def test_hfa_srs_NAD83_CORS96_UTM():
|
|
sr = osr.SpatialReference()
|
|
sr.SetFromUserInput(
|
|
'PROJCS["NAD_1983_CORS96_UTM_Zone_11N",GEOGCS["NAD83(CORS96)",DATUM["NAD83_Continuously_Operating_Reference_Station_1996",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1133"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","6783"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["ESRI","102411"]]'
|
|
)
|
|
|
|
ds = gdal.GetDriverByName("HFA").Create("/vsimem/TestHFASRS.img", 1, 1)
|
|
ds.SetProjection(sr.ExportToWkt())
|
|
ds = None
|
|
|
|
ds = gdal.Open("/vsimem/TestHFASRS.img")
|
|
srs_got = ds.GetSpatialRef()
|
|
assert srs_got.GetAuthorityName(None) is None
|
|
assert srs_got.IsSame(sr), srs_got.ExportToWkt()
|
|
ds = None
|
|
|
|
gdal.Unlink("/vsimem/TestHFASRS.img")
|
|
|
|
|
|
def test_hfa_srs_esri_54049_pe_string_only_broken():
|
|
|
|
with gdaltest.error_handler():
|
|
ds = gdal.Open("../gdrivers/data/hfa/esri_54049_pe_string_only_broken.img")
|
|
assert gdal.GetLastErrorType() == gdal.CE_Warning
|
|
srs_got = ds.GetSpatialRef()
|
|
srs_ref = osr.SpatialReference()
|
|
srs_ref.SetFromUserInput("ESRI:54049")
|
|
assert srs_got.IsSame(srs_ref), srs_got.ExportToWkt()
|
|
|
|
|
|
def test_hfa_srs_DISABLEPESTRING():
|
|
sr = osr.SpatialReference()
|
|
sr.ImportFromEPSG(7844)
|
|
|
|
filename = "/vsimem/test_hfa_srs_DISABLEPESTRING.img"
|
|
ds = gdal.GetDriverByName("HFA").Create(
|
|
filename, 1, 1, options=["DISABLEPESTRING=YES"]
|
|
)
|
|
ds.SetSpatialRef(sr)
|
|
ds = None
|
|
|
|
ds = gdal.Open(filename)
|
|
srs_got = ds.GetSpatialRef()
|
|
# without DISABLEPESTRING, we'd get GCS_GDA2020
|
|
assert srs_got.GetName() == "Geocentric_Datum_of_Australia_2020"
|
|
ds = None
|
|
|
|
gdal.Unlink(filename)
|