Skip to content

Commit

Permalink
ENH: Add support for datum ensemble
Browse files Browse the repository at this point in the history
  • Loading branch information
@snowman2
snowman2 committed Oct 31, 2020
1 parent 340971c commit 23c9d9d
Show file tree
Hide file tree
Showing 6 changed files with 131 additions and 45 deletions.
27 changes: 24 additions & 3 deletions pyproj/_crs.pyx
View file
Expand Up @@ -1288,13 +1288,13 @@ _DATUM_TYPE_MAP = {
}

_PJ_DATUM_TYPE_MAP = {
DatumType.DATUM_ENSEMBLE: PJ_TYPE_DATUM_ENSEMBLE,
DatumType.GEODETIC_REFERENCE_FRAME: PJ_TYPE_GEODETIC_REFERENCE_FRAME,
DatumType.DYNAMIC_GEODETIC_REFERENCE_FRAME:
PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME,
DatumType.VERTICAL_REFERENCE_FRAME: PJ_TYPE_VERTICAL_REFERENCE_FRAME,
DatumType.DYNAMIC_VERTICAL_REFERENCE_FRAME:
PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME,
DatumType.DATUM_ENSEMBLE: PJ_TYPE_DATUM_ENSEMBLE,
}


Expand Down Expand Up @@ -1329,7 +1329,7 @@ cdef class Datum(_CRSParts):
return datum

@staticmethod
def from_authority(auth_name, code):
def _from_authority(auth_name, code, PJ_CATEGORY category):
"""
Create a Datum from an authority code.
Expand All @@ -1350,7 +1350,7 @@ cdef class Datum(_CRSParts):
context,
cstrencode(str(auth_name)),
cstrencode(str(code)),
PJ_CATEGORY_DATUM,
category,
False,
NULL,
)
Expand All @@ -1361,6 +1361,27 @@ cdef class Datum(_CRSParts):
CRSError.clear()
return Datum.create(context, datum_pj)

@staticmethod
def from_authority(auth_name, code):
"""
Create a Datum from an authority code.
Parameters
----------
auth_name: str
Name ot the authority.
code: str or int
The code used by the authority.
Returns
-------
Datum
"""
try:
return Datum._from_authority(auth_name, code, PJ_CATEGORY_DATUM_ENSEMBLE)
except CRSError:
return Datum._from_authority(auth_name, code, PJ_CATEGORY_DATUM)

@staticmethod
def from_epsg(code):
"""
Expand Down
31 changes: 22 additions & 9 deletions pyproj/proj.pxi
View file
@@ -1,4 +1,25 @@
# PROJ.4 API Defnition
# PROJ API Definition

IF CTE_PROJ_VERSION_MAJOR >= 8:
cdef extern from "proj.h":
ctypedef enum PJ_CATEGORY:
PJ_CATEGORY_ELLIPSOID
PJ_CATEGORY_PRIME_MERIDIAN
PJ_CATEGORY_DATUM
PJ_CATEGORY_CRS
PJ_CATEGORY_COORDINATE_OPERATION
PJ_CATEGORY_DATUM_ENSEMBLE
ELSE:
cdef extern from "proj.h":
ctypedef enum PJ_CATEGORY:
PJ_CATEGORY_ELLIPSOID
PJ_CATEGORY_PRIME_MERIDIAN
PJ_CATEGORY_DATUM
PJ_CATEGORY_CRS
PJ_CATEGORY_COORDINATE_OPERATION
cdef int PJ_CATEGORY_DATUM_ENSEMBLE = PJ_CATEGORY_DATUM


cdef extern from "proj.h":
cdef int PROJ_VERSION_MAJOR
cdef int PROJ_VERSION_MINOR
Expand Down Expand Up @@ -397,20 +418,12 @@ cdef extern from "proj.h":
PJ *proj_concatoperation_get_step(PJ_CONTEXT *ctx,
const PJ *concatoperation,
int i_step)

ctypedef enum PJ_CATEGORY:
PJ_CATEGORY_ELLIPSOID
PJ_CATEGORY_PRIME_MERIDIAN
PJ_CATEGORY_DATUM
PJ_CATEGORY_CRS
PJ_CATEGORY_COORDINATE_OPERATION
PJ *proj_create_from_database(PJ_CONTEXT *ctx,
const char *auth_name,
const char *code,
PJ_CATEGORY category,
int usePROJAlternativeGridNames,
const char* const *options)

PJ_OBJ_LIST *proj_create_from_name(PJ_CONTEXT *ctx,
const char *auth_name,
const char *searchedName,
Expand Down
35 changes: 22 additions & 13 deletions setup.py
View file
Expand Up @@ -14,23 +14,26 @@
INTERNAL_PROJ_DIR = CURRENT_FILE_PATH / "pyproj" / BASE_INTERNAL_PROJ_DIR


def check_proj_version(proj_dir: Path):
"""checks that the PROJ library meets the minimum version"""
def get_proj_version(proj_dir: Path) -> str:
proj_version = os.environ.get("PROJ_VERSION")
if proj_version:
return proj_version
proj = proj_dir / "bin" / "proj"
proj_ver_bytes = subprocess.check_output(
str(proj), stderr=subprocess.STDOUT
).decode("ascii")
proj_ver_bytes = (proj_ver_bytes.split()[1]).strip(",")
proj_version = parse_version(proj_ver_bytes)
if proj_version < PROJ_MIN_VERSION:
proj_ver = subprocess.check_output(str(proj), stderr=subprocess.STDOUT).decode(
"ascii"
)
return (proj_ver.split()[1]).strip(",")


def check_proj_version(proj_version: str) -> None:
"""checks that the PROJ library meets the minimum version"""
if parse_version(proj_version) < PROJ_MIN_VERSION:
raise SystemExit(
f"ERROR: Minimum supported proj version is {PROJ_MIN_VERSION}, installed "
f"version is {proj_version}. For more information see: "
"https://pyproj4.github.io/pyproj/stable/installation.html"
)

return proj_version


def get_proj_dir() -> Path:
"""
Expand Down Expand Up @@ -58,9 +61,6 @@ def get_proj_dir() -> Path:
print("PROJ_DIR is set, using existing proj4 installation..\n")
else:
raise SystemExit(f"ERROR: Invalid path for PROJ_DIR {proj_dir}")

# check_proj_version
check_proj_version(proj_dir)
return proj_dir


Expand Down Expand Up @@ -152,6 +152,10 @@ def get_extension_modules():
library_dirs = get_proj_libdirs(proj_dir)
include_dirs = get_proj_incdirs(proj_dir)

proj_version = get_proj_version(proj_dir)
check_proj_version(proj_version)
proj_version_major, proj_version_minor, proj_version_patch = proj_version.split(".")

# setup extension options
ext_options = {
"include_dirs": include_dirs,
Expand All @@ -175,6 +179,11 @@ def get_extension_modules():
Extension("pyproj._sync", ["pyproj/_sync.pyx"], **ext_options),
],
quiet=True,
compile_time_env={
"CTE_PROJ_VERSION_MAJOR": int(proj_version_major),
"CTE_PROJ_VERSION_MINOR": int(proj_version_minor),
"CTE_PROJ_VERSION_PATCH": int(proj_version_patch),
},
**get_cythonize_options(),
)

Expand Down
8 changes: 8 additions & 0 deletions test/conftest.py
View file
@@ -1,11 +1,13 @@
import os
from contextlib import contextmanager
from distutils.version import LooseVersion
from pathlib import Path

import pyproj
from pyproj.datadir import get_data_dir, get_user_data_dir, set_data_dir

_NETWORK_ENABLED = pyproj.network.is_network_enabled()
PROJ_GTE_8 = LooseVersion(pyproj.__proj_version__) >= LooseVersion("8.0")


def unset_data_dir():
Expand Down Expand Up @@ -70,3 +72,9 @@ def grids_available(*grid_names, check_network=True, check_all=False):
if check_all:
return all(available)
return any(available)


def get_wgs84_datum_name():
if PROJ_GTE_8:
return "World Geodetic System 1984 ensemble"
return "World Geodetic System 1984"
72 changes: 53 additions & 19 deletions test/crs/test_crs.py
View file
Expand Up @@ -15,7 +15,7 @@
from pyproj.enums import ProjVersion, WktVersion
from pyproj.exceptions import CRSError
from pyproj.transformer import TransformerGroup
from test.conftest import grids_available
from test.conftest import PROJ_GTE_8, get_wgs84_datum_name, grids_available


class CustomCRS(object):
Expand Down Expand Up @@ -201,25 +201,28 @@ def test_repr():
"Area of Use:\n"
"- name: World.\n"
"- bounds: (-180.0, -90.0, 180.0, 90.0)\n"
"Datum: World Geodetic System 1984\n"
f"Datum: {get_wgs84_datum_name()}\n"
"- Ellipsoid: WGS 84\n"
"- Prime Meridian: Greenwich\n"
)


def test_repr__long():
with pytest.warns(FutureWarning):
if PROJ_GTE_8:
wkt_str = 'GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1'
else:
wkt_str = 'GEOGCRS["WGS 84",DATUM["World Geodetic System 1984'
assert repr(CRS(CRS({"init": "EPSG:4326"}).to_wkt())) == (
'<Geographic 2D CRS: GEOGCRS["WGS 84",'
'DATUM["World Geodetic System 1984 ...>\n'
f"<Geographic 2D CRS: {wkt_str} ...>\n"
"Name: WGS 84\n"
"Axis Info [ellipsoidal]:\n"
"- lon[east]: Longitude (degree)\n"
"- lat[north]: Latitude (degree)\n"
"Area of Use:\n"
"- name: World.\n"
"- bounds: (-180.0, -90.0, 180.0, 90.0)\n"
"Datum: World Geodetic System 1984\n"
f"Datum: {get_wgs84_datum_name()}\n"
"- Ellipsoid: WGS 84\n"
"- Prime Meridian: Greenwich\n"
)
Expand All @@ -235,7 +238,7 @@ def test_repr_epsg():
"Area of Use:\n"
"- name: World.\n"
"- bounds: (-180.0, -90.0, 180.0, 90.0)\n"
"Datum: World Geodetic System 1984\n"
f"Datum: {get_wgs84_datum_name()}\n"
"- Ellipsoid: WGS 84\n"
"- Prime Meridian: Greenwich\n"
)
Expand Down Expand Up @@ -309,9 +312,13 @@ def test_epsg():

def test_datum():
datum = CRS.from_epsg(4326).datum
assert repr(datum).startswith('DATUM["World Geodetic System 1984"')
assert "\n" in repr(datum)
assert datum.to_wkt().startswith('DATUM["World Geodetic System 1984"')
if PROJ_GTE_8:
datum_wkt = 'ENSEMBLE["World Geodetic System 1984 ensemble"'
else:
datum_wkt = 'DATUM["World Geodetic System 1984"'
assert repr(datum).startswith(datum_wkt)
assert datum.to_wkt().startswith(datum_wkt)
assert datum == datum
assert datum.is_exact_same(datum)

Expand Down Expand Up @@ -621,16 +628,31 @@ def test_coordinate_operation__from_authority__empty():


def test_datum__from_epsg():
assert Datum.from_epsg("6326").to_wkt() == (
'DATUM["World Geodetic System 1984",'
'ELLIPSOID["WGS 84",6378137,298.257223563,'
'LENGTHUNIT["metre",1]],ID["EPSG",6326]]'
)
if PROJ_GTE_8:
datum_wkt = (
'ENSEMBLE["World Geodetic System 1984 ensemble",'
'MEMBER["World Geodetic System 1984 (Transit)",'
'ID["EPSG",1166]],MEMBER["World Geodetic System 1984 (G730)",'
'ID["EPSG",1152]],MEMBER["World Geodetic System 1984 (G873)",'
'ID["EPSG",1153]],MEMBER["World Geodetic System 1984 (G1150)",'
'ID["EPSG",1154]],MEMBER["World Geodetic System 1984 (G1674)",'
'ID["EPSG",1155]],MEMBER["World Geodetic System 1984 (G1762)",'
'ID["EPSG",1156]],ELLIPSOID["WGS 84",6378137,298.257223563,'
'LENGTHUNIT["metre",1],ID["EPSG",7030]],'
'ENSEMBLEACCURACY[2.0],ID["EPSG",6326]]'
)
else:
datum_wkt = (
'DATUM["World Geodetic System 1984",'
'ELLIPSOID["WGS 84",6378137,298.257223563,'
'LENGTHUNIT["metre",1]],ID["EPSG",6326]]'
)
assert Datum.from_epsg("6326").to_wkt() == datum_wkt


def test_datum__from_authority():
dt = Datum.from_authority("EPSG", 6326)
assert dt.name == "World Geodetic System 1984"
assert dt.name == get_wgs84_datum_name()


def test_datum__from_epsg__invalid():
Expand All @@ -648,7 +670,9 @@ def test_datum__from_authority__invalid():
[
6326,
("EPSG", 6326),
"urn:ogc:def:datum:EPSG::6326",
"urn:ogc:def:ensemble:EPSG::6326"
if PROJ_GTE_8
else "urn:ogc:def:datum:EPSG::6326",
Datum.from_epsg(6326),
Datum.from_epsg(6326).to_json_dict(),
"World Geodetic System 1984",
Expand Down Expand Up @@ -867,12 +891,22 @@ def test_datum_equals():


@pytest.mark.parametrize(
"input_str", ["urn:ogc:def:datum:EPSG::6326", "World Geodetic System 1984"]
"input_str",
[
"urn:ogc:def:ensemble:EPSG::6326"
if PROJ_GTE_8
else "urn:ogc:def:datum:EPSG::6326",
"World Geodetic System 1984",
],
)
def test_datum__from_string(input_str):
dd = Datum.from_string(input_str)
assert dd.name == "World Geodetic System 1984"
assert dd.type_name == "Geodetic Reference Frame"
if PROJ_GTE_8:
assert dd.name == "World Geodetic System 1984 ensemble"
assert dd.type_name == "Datum Ensemble"
else:
assert dd.name == "World Geodetic System 1984"
assert dd.type_name == "Geodetic Reference Frame"


@pytest.mark.parametrize(
Expand All @@ -897,7 +931,7 @@ def test_datum__from_string__type_name(input_str, type_name):
)
def test_datum__from_name(input_name):
dd = Datum.from_name(input_name)
assert dd.name == "World Geodetic System 1984"
assert dd.name == get_wgs84_datum_name()


@pytest.mark.parametrize("auth_name", [None, "ESRI"])
Expand Down
3 changes: 2 additions & 1 deletion test/crs/test_crs_cf.py
View file
Expand Up @@ -16,6 +16,7 @@
VerticalPerspectiveConversion,
)
from pyproj.exceptions import CRSError
from test.conftest import get_wgs84_datum_name


def _to_dict(operation):
Expand Down Expand Up @@ -240,7 +241,7 @@ def test_cf_from_utm():
"longitude_of_prime_meridian": 0.0,
"prime_meridian_name": "Greenwich",
"geographic_crs_name": "WGS 84",
"horizontal_datum_name": "World Geodetic System 1984",
"horizontal_datum_name": get_wgs84_datum_name(),
"projected_crs_name": "WGS 84 / UTM zone 15N",
"grid_mapping_name": "transverse_mercator",
"latitude_of_projection_origin": 0.0,
Expand Down

0 comments on commit 23c9d9d

Please sign in to comment.