Authors: SeanGillies (SG), HowardButler? (HB)
A crucial part of the library is a class representing a spatial reference system?.
There are three primary ways of defining SRS:
1. So-called WellKnownText?, such as:
PROJCS["NAD83(HARN) / Colorado North (ftUS)",
GEOGCS["NAD83(HARN)",
DATUM["NAD83_High_Accuracy_Regional_Network",
SPHEROID["GRS 1980", 6378137, 298.257222101,
AUTHORITY["EPSG", "7019"]],
AUTHORITY["EPSG", "6152"]],
PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]],
AUTHORITY["EPSG", "4152"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1", 40.78333333333333],
PARAMETER["standard_parallel_2", 39.71666666666667],
PARAMETER["latitude_of_origin", 39.33333333333334],
PARAMETER["central_meridian", -105.5],
PARAMETER["false_easting", 3000000],
PARAMETER["false_northing", 1000000],
UNIT["US survey foot", 0.3048006096012192, AUTHORITY["EPSG", "9003"]],
AUTHORITY["EPSG", "2876"]]
--There are at least three variants of WKT. ESRI, OGC, and Oracle. ESRI WKT is *required* for sucking in .prj file coordinate systems. (HB)
2. EPSG codes: 4326 being the long/lat WGS84, and 26913 being NAD83 UTM zone 13. 2876 is the code for the projection above. Note how that appears in the top-level authority.
3. PROJ.4 parameterization:
+proj=lcc +lat_1=40.78333333333333 +lat_2=39.71666666666667 \ +lat_0=39.33333333333334 +lon_0=-105.5 +x_0=914401.8288036576 \ +y_0=304800.6096012192 +ellps=GRS80 +to_meter=0.3048006096012192 \ +no_defs
PCL depends on the PROJ.4 library for transformations, and so it makes sense to keep close to that parameterization as our model. Conversion to and from the other forms is desirable.
After a chat with Howard, I think the default should be initialization from a proj4 string:
srs = SpatialReference('+proj=lcc +lat_1=...')
and that we shall keep the EPSG shortcut
srs = SpatialReference(epsg=2876)
and add a WKT option
srs = SpatialReference(wkt='PROJCS[...]')
We could also initialize from proj4-style keyword args:
srs = SpatialReference(proj='lcc', lat_1=40.783, ...)
--Python's indefinite precision, ie a 1.10 is actually 1.10000000000001, could cause significant problems if we initialize this way. (HB).
SG: do you mean problems with comparing for equality with a projection in the epsg table?
This would suggest an initializer method like
class SpatialReference:
def __init__(self, proj4=None, epsg=None, wkt=None, **kw):
...
Clarified to answer a question by Howard.
Instances of SpatialReference? should probably be, in practice, immutable after creation. I don't see the need to modify an SRS. Simply discard it and create a new one.
SRS Translation
IMO, it is very important to be able to translate one SRS format to another. GDAL/OGR provides this via the OSR API. This allows us to take in WKT in ESRI format and export it as OGC WKT, or take in an EPSG code and emit ESRI formatted SRS. Here's a C Python module that uses OSR to do this (HB):
/*
# =============================================================================
# OSR Reference Utilities Module Copyright(C) 2006 Howard C. Butler
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 2 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 59 Temple
# Place, Suite 330, Boston, MA 02111-1307 USA
#
# Contact email: hobu@hobu.net
# =============================================================================
*/
#include "Python.h"
#include <ogr/ogr_srs_api.h>
#include <port/cpl_port.h>
#include <port/cpl_error.h>
static PyObject *OSRError;
void ThrowIfCPLError(OGRErr err)
{
if (err)
{
PyErr_SetString(OSRError, CPLGetLastErrorMsg());
CPLErrorReset();
}
}
PyObject* to_wkt(OGRSpatialReferenceH *ref)
{
char* wkt;
OGRErr err;
PyObject* out_string;
err = OSRExportToWkt(ref, &wkt);
if (!err)
{
out_string= PyString_FromString(wkt);
free(wkt);
return out_string;
}
else
{
ThrowIfCPLError(err);
return NULL;
}
}
PyObject* to_xml(OGRSpatialReferenceH *ref)
{
char* wkt;
OGRErr err;
PyObject* out_string;
err = OSRExportToXML(ref, &wkt, NULL);
if (!err)
{
out_string= PyString_FromString(wkt);
free(wkt);
return out_string;
}
else
{
ThrowIfCPLError(err);
return NULL;
}
}
PyObject* to_proj4(OGRSpatialReferenceH *ref)
{
char* wkt;
OGRErr err;
PyObject* out_string;
err = OSRExportToProj4(ref, &wkt);
if (!err)
{
out_string= PyString_FromString(wkt);
free(wkt);
return out_string;
}
else
{
ThrowIfCPLError(err);
return NULL;
}
}
PyObject* to_esri(OGRSpatialReferenceH *ref)
{
char* wkt;
OGRErr err;
PyObject* out_string;
err = OSRMorphToESRI(ref);
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
err = OSRExportToWkt(ref, &wkt);
if (!err)
{
out_string= PyString_FromString(wkt);
free(wkt);
return out_string;
}
else
{
ThrowIfCPLError(err);
return NULL;
}
}
OGRSpatialReferenceH * from_esri(char* ref_defn)
{
OGRErr err;
OGRSpatialReferenceH ref=NULL;
char *list[2] = {NULL,NULL};
ref = OSRNewSpatialReference(NULL);
list[0] = ref_defn;
list[1] = NULL;
err = OSRImportFromESRI(ref, list);
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
return ref;
}
OGRSpatialReferenceH * from_epsg(char* ref_defn)
{
OGRErr err;
OGRSpatialReferenceH ref=NULL;
ref = OSRNewSpatialReference(NULL);
err = OSRImportFromEPSG(ref, atoi(ref_defn));
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
return ref;
}
OGRSpatialReferenceH * from_proj4(char* ref_defn)
{
OGRErr err;
OGRSpatialReferenceH ref=NULL;
ref = OSRNewSpatialReference(NULL);
err = OSRImportFromProj4(ref, ref_defn);
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
return ref;
}
OGRSpatialReferenceH * from_xml(char* ref_defn)
{
OGRErr err;
OGRSpatialReferenceH ref=NULL;
ref = OSRNewSpatialReference(NULL);
err = OSRImportFromXML(ref, ref_defn);
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
return ref;
}
OGRSpatialReferenceH * from_wkt(char* ref_defn)
{
OGRErr err;
OGRSpatialReferenceH ref=NULL;
ref = OSRNewSpatialReference(NULL);
err = OSRImportFromWkt(ref, &ref_defn);
if (err)
{
ThrowIfCPLError(err);
return NULL;
}
return ref;
}
int validate_input(char* input)
{
if (!strcmp(input,"epsg"))
return 1;
if (!strcmp(input,"esri"))
return 1;
if (!strcmp(input,"proj4"))
return 1;
if (!strcmp(input,"wkt"))
return 1;
if (!strcmp(input,"xml"))
return 1;
return 0;
}
int validate_output(char* input)
{
if (!strcmp(input,"esri"))
return 1;
if (!strcmp(input,"proj4"))
return 1;
if (!strcmp(input,"wkt"))
return 1;
if (!strcmp(input,"xml"))
return 1;
return 0;
}
static PyObject *
_translate(PyObject *self, PyObject *args)
{
PyObject *output=NULL;
char* in_ref_type;
char* out_ref_type;
char* ref_defn;
OGRSpatialReferenceH ref=NULL;
if (!PyArg_ParseTuple(args, "sss:_translate", &ref_defn, &in_ref_type, &out_ref_type))
goto fail;
if (!validate_input(in_ref_type))
{
PyErr_SetString(OSRError, "input SRS type not recognized ... use 'esri','proj4','wkt','xml', or 'epsg'");
goto fail;
}
if (!validate_output(out_ref_type))
{
PyErr_SetString(OSRError, "output SRS type not recognized... use 'esri','proj4','wkt', or 'xml'");
goto fail;
}
if (!strcmp(in_ref_type,"epsg"))
{
ref = from_epsg(ref_defn);
}
if (!strcmp(in_ref_type,"esri"))
{
ref = from_esri(ref_defn);
}
if (!strcmp(in_ref_type,"proj4"))
{
ref = from_proj4(ref_defn);
}
if (!strcmp(in_ref_type,"xml"))
{
ref = from_xml(ref_defn);
}
if (!strcmp(in_ref_type,"wkt"))
{
ref = from_wkt(ref_defn);
}
if (ref)
{
if (!strcmp(out_ref_type,"wkt"))
{
output = to_wkt(ref);
}
if (!strcmp(out_ref_type,"xml"))
{
output = to_xml(ref);
}
if (!strcmp(out_ref_type,"esri"))
{
output = to_esri(ref);
}
if (!strcmp(out_ref_type,"proj4"))
{
output = to_proj4(ref);
}
}
if (output)
{
return output;
}
else
{
return NULL;
}
fail:
return NULL;
}
/*
===============================================================================
Initialization
===============================================================================
*/
static PyMethodDef OSRMethods[] = {
{"_translate", (PyCFunction)_translate, METH_VARARGS,
"Translates an SRS definition in one format to another"},
{NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC
init_osrtranslate(void)
{
PyObject* module;
module = Py_InitModule3("_osrtranslate", OSRMethods,
"OSR reference system translation module");
if (module == NULL)
return;
OSRError = PyErr_NewException("_osrtranslate.OSRError", NULL, NULL);
Py_INCREF(OSRError);
PyModule_AddObject(module, "OSRError", OSRError);
/* Suppress GDAL from emitting errors to stdout by default */
CPLPushErrorHandler(CPLQuietErrorHandler);
}
