Loading...
Searching...
No Matches
convert.c
Go to the documentation of this file.
/*!
\file lib/proj/convert.c
\brief GProj Library - Functions for manipulating co-ordinate
system representations
(C) 2003-2018 by the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
\author Paul Kelly, Frank Warmerdam, Markus Metz
*/
#include <grass/config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
#include <cpl_csv.h>
#include "local_proto.h"
/* GRASS relative location of OGR co-ordinate system lookup tables */
#define CSVDIR "/etc/proj/ogr_csv"
static void DatumNameMassage(char **);
/* from proj-5.0.0/src/pj_units.c */
struct gpj_units {
char *id; /* units keyword */
char *to_meter; /* multiply by value to get meters */
char *name; /* comments */
double factor; /* to_meter factor in actual numbers */
};
struct gpj_units gpj_units[] = {
{"km", "1000.", "Kilometer", 1000.0},
{"m", "1.", "Meter", 1.0},
{"dm", "1/10", "Decimeter", 0.1},
{"cm", "1/100", "Centimeter", 0.01},
{"mm", "1/1000", "Millimeter", 0.001},
{"kmi", "1852.0", "International Nautical Mile", 1852.0},
{"in", "0.0254", "International Inch", 0.0254},
{"ft", "0.3048", "International Foot", 0.3048},
{"yd", "0.9144", "International Yard", 0.9144},
{"mi", "1609.344", "International Statute Mile", 1609.344},
{"fath", "1.8288", "International Fathom", 1.8288},
{"ch", "20.1168", "International Chain", 20.1168},
{"link", "0.201168", "International Link", 0.201168},
{"us-in", "1./39.37", "U.S. Surveyor's Inch", 0.0254},
{"us-ft", "0.304800609601219", "U.S. Surveyor's Foot", 0.304800609601219},
{"us-yd", "0.914401828803658", "U.S. Surveyor's Yard", 0.914401828803658},
{"us-ch", "20.11684023368047", "U.S. Surveyor's Chain", 20.11684023368047},
{"us-mi", "1609.347218694437", "U.S. Surveyor's Statute Mile",
1609.347218694437},
{"ind-yd", "0.91439523", "Indian Yard", 0.91439523},
{"ind-ft", "0.30479841", "Indian Foot", 0.30479841},
{"ind-ch", "20.11669506", "Indian Chain", 20.11669506},
{NULL, NULL, NULL, 0.0}};
static char *grass_to_wkt(const struct Key_Value *proj_info,
const struct Key_Value *proj_units,
const struct Key_Value *proj_epsg, int esri_style,
int prettify)
{
OGRSpatialReferenceH hSRS;
char *wkt, *local_wkt;
hSRS = GPJ_grass_to_osr2(proj_info, proj_units, proj_epsg);
if (hSRS == NULL)
return NULL;
if (esri_style)
OSRMorphToESRI(hSRS);
if (prettify)
OSRExportToPrettyWkt(hSRS, &wkt, 0);
else
OSRExportToWkt(hSRS, &wkt);
local_wkt = G_store(wkt);
CPLFree(wkt);
OSRDestroySpatialReference(hSRS);
return local_wkt;
}
/*!
* \brief Converts a GRASS co-ordinate system representation to WKT style.
*
* Takes a GRASS co-ordinate system as specified by two sets of
* key/value pairs derived from the PROJ_INFO and PROJ_UNITS files,
* and converts it to the 'Well Known Text' format.
*
* \param proj_info Set of GRASS PROJ_INFO key/value pairs
* \param proj_units Set of GRASS PROJ_UNIT key/value pairs
* \param esri_style boolean Output ESRI-style WKT (Use OSRMorphToESRI()
* function provided by OGR library)
* \param prettify boolean Use linebreaks and indents to 'prettify' output
* WKT string (Use OSRExportToPrettyWkt() function in OGR)
*
* \return Pointer to a string containing the co-ordinate system in
* WKT format
* \return NULL on error
*/
char *GPJ_grass_to_wkt(const struct Key_Value *proj_info,
const struct Key_Value *proj_units, int esri_style,
int prettify)
{
return grass_to_wkt(proj_info, proj_units, NULL, esri_style, prettify);
}
/*!
* \brief Converts a GRASS co-ordinate system representation to WKT
* style. EPSG code is preferred if available.
*
* Takes a GRASS co-ordinate system as specified key/value pairs
* derived from the PROJ_EPSG file. TOWGS84 parameter is scanned
* from PROJ_INFO file and appended to co-ordinate system definition
* imported from EPSG code by GDAL library. PROJ_UNITS file is
* ignored. The function converts it to the 'Well Known Text' format.
*
* \todo Merge with GPJ_grass_to_wkt() in GRASS 8.
*
* \param proj_info Set of GRASS PROJ_INFO key/value pairs
* \param proj_units Set of GRASS PROJ_UNIT key/value pairs
* \param proj_epsg Set of GRASS PROJ_EPSG key/value pairs
* \param esri_style boolean Output ESRI-style WKT (Use OSRMorphToESRI()
* function provided by OGR library)
* \param prettify boolean Use linebreaks and indents to 'prettify' output
* WKT string (Use OSRExportToPrettyWkt() function in OGR)
*
* \return Pointer to a string containing the co-ordinate system in
* WKT format
* \return NULL on error
*/
char *GPJ_grass_to_wkt2(const struct Key_Value *proj_info,
const struct Key_Value *proj_units,
const struct Key_Value *proj_epsg, int esri_style,
int prettify)
{
return grass_to_wkt(proj_info, proj_units, proj_epsg, esri_style, prettify);
}
/*!
* \brief Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
*
* \param proj_info Set of GRASS PROJ_INFO key/value pairs
* \param proj_units Set of GRASS PROJ_UNIT key/value pairs
*
* \return OGRSpatialReferenceH object representing the co-ordinate system
* defined by proj_info and proj_units or NULL if it fails
*/
OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value *proj_info,
const struct Key_Value *proj_units)
{
struct pj_info pjinfo;
char *proj4, *proj4mod, *wkt, *modwkt, *startmod, *lastpart;
OGRSpatialReferenceH hSRS, hSRS2;
OGRErr errcode;
struct gpj_datum dstruct;
struct gpj_ellps estruct;
size_t len;
const char *ellpskv, *unit, *unfact;
char *ellps, *ellpslong, *datum, *params, *towgs84, *datumlongname, *start,
*end;
const char *sysname, *osrunit;
double a, es, rf;
int haveparams = 0;
if ((proj_info == NULL) || (proj_units == NULL))
return NULL;
hSRS = OSRNewSpatialReference(NULL);
/* create PROJ structure from GRASS key/value pairs */
if (pj_get_kv(&pjinfo, proj_info, proj_units) < 0) {
G_warning(_("Unable parse GRASS PROJ_INFO file"));
return NULL;
}
/* fetch the PROJ definition */
/* TODO: get the PROJ definition as used by pj_get_kv() */
if ((proj4 = pjinfo.def) == NULL) {
G_warning(_("Unable get PROJ.4-style parameter string"));
return NULL;
}
proj_destroy(pjinfo.pj);
unit = G_find_key_value("unit", proj_units);
unfact = G_find_key_value("meters", proj_units);
if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
else
proj4mod = G_store(proj4);
/* create GDAL OSR from proj string */
if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
G_warning(_("OGR can't parse PROJ.4-style parameter string: "
"%s (OGR Error code was %d)"),
proj4mod, errcode);
return NULL;
}
G_free(proj4mod);
/* this messes up PROJCS versus GEOGCS!
sysname = G_find_key_value("name", proj_info);
if (sysname)
OSRSetProjCS(hSRS, sysname);
*/
if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
G_warning(_("OGR can't get WKT-style parameter string "
"(OGR Error code was %d)"),
errcode);
return NULL;
}
ellpskv = G_find_key_value("ellps", proj_info);
GPJ__get_ellipsoid_params(proj_info, &a, &es, &rf);
haveparams = GPJ__get_datum_params(proj_info, &datum, &params);
if (ellpskv != NULL)
ellps = G_store(ellpskv);
else
ellps = NULL;
if ((datum == NULL) || (GPJ_get_datum_by_name(datum, &dstruct) < 0)) {
datumlongname = G_store("unknown");
if (ellps == NULL)
ellps = G_store("unnamed");
}
else {
datumlongname = G_store(dstruct.longname);
if (ellps == NULL)
ellps = G_store(dstruct.ellps);
GPJ_free_datum(&dstruct);
}
G_debug(3, "GPJ_grass_to_osr: datum: <%s>", datum);
G_free(datum);
if (GPJ_get_ellipsoid_by_name(ellps, &estruct) > 0) {
ellpslong = G_store(estruct.longname);
DatumNameMassage(&ellpslong);
GPJ_free_ellps(&estruct);
}
else
ellpslong = G_store(ellps);
startmod = strstr(wkt, "GEOGCS");
lastpart = strstr(wkt, "PRIMEM");
len = strlen(wkt) - strlen(startmod);
wkt[len] = '\0';
if (haveparams == 2) {
/* Only put datum params into the WKT if they were specifically
* specified in PROJ_INFO */
char *paramkey, *paramvalue;
paramkey = strtok(params, "=");
paramvalue = params + strlen(paramkey) + 1;
if (G_strcasecmp(paramkey, "towgs84") == 0)
G_asprintf(&towgs84, ",TOWGS84[%s]", paramvalue);
else
towgs84 = G_store("");
G_free(params);
}
else
towgs84 = G_store("");
sysname = OSRGetAttrValue(hSRS, "PROJCS", 0);
if (sysname == NULL) {
/* Not a projected co-ordinate system */
start = G_store("");
end = G_store("");
}
else {
if ((strcmp(sysname, "unnamed") == 0) &&
(G_find_key_value("name", proj_info) != NULL))
G_asprintf(&start, "PROJCS[\"%s\",",
G_find_key_value("name", proj_info));
else
start = G_store(wkt);
osrunit = OSRGetAttrValue(hSRS, "UNIT", 0);
if ((unfact == NULL) || (G_strcasecmp(osrunit, "unknown") != 0))
end = G_store("");
else {
char *buff;
double unfactf = atof(unfact);
G_asprintf(&buff, ",UNIT[\"%s\",", osrunit);
startmod = strstr(lastpart, buff);
len = strlen(lastpart) - strlen(startmod);
lastpart[len] = '\0';
G_free(buff);
if (unit == NULL)
unit = "unknown";
G_asprintf(&end, ",UNIT[\"%s\",%.16g]]", unit, unfactf);
}
}
OSRDestroySpatialReference(hSRS);
G_asprintf(
&modwkt,
"%sGEOGCS[\"%s\",DATUM[\"%s\",SPHEROID[\"%s\",%.16g,%.16g]%s],%s%s",
start, ellps, datumlongname, ellpslong, a, rf, towgs84, lastpart, end);
hSRS2 = OSRNewSpatialReference(modwkt);
G_free(modwkt);
CPLFree(wkt);
G_free(start);
G_free(ellps);
G_free(datumlongname);
G_free(ellpslong);
G_free(towgs84);
G_free(end);
return hSRS2;
}
/*!
* \brief Converts a GRASS co-ordinate system to an
* OGRSpatialReferenceH object. EPSG code is preferred if available.
*
* The co-ordinate system definition is imported from EPSG (by GDAL)
* definition if available. TOWGS84 parameter is scanned from
* PROJ_INFO file and appended to co-ordinate system definition. If
* EPSG code is not available, PROJ_INFO file is used as
* GPJ_grass_to_osr() does.
* \todo Merge with GPJ_grass_to_osr() in GRASS 8.
*
* \param proj_info Set of GRASS PROJ_INFO key/value pairs
* \param proj_units Set of GRASS PROJ_UNIT key/value pairs
* \param proj_epsg Set of GRASS PROJ_EPSG key/value pairs
*
* \return OGRSpatialReferenceH object representing the co-ordinate system
* defined by proj_info and proj_units or NULL if it fails
*/
OGRSpatialReferenceH GPJ_grass_to_osr2(const struct Key_Value *proj_info,
const struct Key_Value *proj_units,
const struct Key_Value *proj_epsg)
{
int epsgcode = 0;
if (proj_epsg) {
const char *epsgstr = G_find_key_value("epsg", proj_epsg);
if (epsgstr)
epsgcode = atoi(epsgstr);
}
if (epsgcode) {
const char *towgs84;
OGRSpatialReferenceH hSRS;
hSRS = OSRNewSpatialReference(NULL);
OSRImportFromEPSG(hSRS, epsgcode);
/* take +towgs84 from projinfo file if defined */
towgs84 = G_find_key_value("towgs84", proj_info);
if (towgs84) {
char **tokens;
int i;
double df[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
tokens = G_tokenize(towgs84, ",");
for (i = 0; i < G_number_of_tokens(tokens); i++)
df[i] = atof(tokens[i]);
G_free_tokens(tokens);
OSRSetTOWGS84(hSRS, df[0], df[1], df[2], df[3], df[4], df[5],
df[6]);
}
return hSRS;
}
return GPJ_grass_to_osr(proj_info, proj_units);
}
/*!
* \brief Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
*
* \param cellhd Pointer to a GRASS Cell_head structure that will have its
* projection-related members populated with appropriate
* values
* \param projinfo Pointer to a pointer which will have a GRASS
* Key_Value structure allocated containing a set of GRASS
* PROJ_INFO values
* \param projunits Pointer to a pointer which will have a GRASS Key_Value
* structure allocated containing a set of GRASS PROJ_UNITS
* values.
* \param hSRS1 OGRSpatialReferenceH object containing the
* co-ordinate system to be converted
* \param datumtrans Index number of datum parameter set to use, 0 to leave
* unspecified
*
* \return 2 if a projected or lat/long co-ordinate system has been
* defined
* \return 1 if an unreferenced XY co-ordinate system has
* been defined
*/
int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
struct Key_Value **projunits, OGRSpatialReferenceH hSRS1,
int datumtrans)
{
struct Key_Value *temp_projinfo, *temp_projinfo_ext;
char *pszProj4 = NULL, *pszRemaining;
char *pszProj = NULL;
const char *pszProjCS = NULL;
char *datum = NULL;
char *proj4_unit = NULL;
struct gpj_datum dstruct;
const char *ograttr;
OGRSpatialReferenceH hSRS;
int use_proj_extension;
*projinfo = NULL;
*projunits = NULL;
hSRS = hSRS1;
if (hSRS == NULL)
goto default_to_xy;
/* Set finder function for locating OGR csv co-ordinate system tables */
/* SetCSVFilenameHook(GPJ_set_csv_loc); */
/* Hopefully this doesn't do any harm if it wasn't in ESRI format
* to start with... */
OSRMorphFromESRI(hSRS);
*projinfo = G_create_key_value();
use_proj_extension = 0;
/* use proj4 definition from EXTENSION attribute if existing */
ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 0);
if (ograttr && *ograttr && strcmp(ograttr, "PROJ4") == 0) {
ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 1);
G_debug(3, "proj4 extension:");
G_debug(3, "%s", ograttr);
if (ograttr && *ograttr) {
char *proj4ext;
OGRSpatialReferenceH hSRS2;
hSRS2 = OSRNewSpatialReference(NULL);
proj4ext = G_store(ograttr);
/* test */
if (OSRImportFromProj4(hSRS2, proj4ext) != OGRERR_NONE) {
G_warning(_("Updating spatial reference with embedded proj4 "
"definition failed. "
"Proj4 definition: <%s>"),
proj4ext);
OSRDestroySpatialReference(hSRS2);
}
else {
/* use new OGR spatial reference defined with embedded proj4
* string */
/* TODO: replace warning with important_message once confirmed
* working */
G_warning(_("Updating spatial reference with embedded proj4 "
"definition"));
/* --------------------------------------------------------------------
*/
/* Derive the user name for the coordinate system. */
/* --------------------------------------------------------------------
*/
pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
if (!pszProjCS)
pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
if (pszProjCS) {
G_set_key_value("name", pszProjCS, *projinfo);
}
else if (pszProj) {
char path[4095];
char name[80];
/* use name of the projection as name for the coordinate
* system */
snprintf(path, sizeof(path), "%s/etc/proj/projections",
G_gisbase());
if (G_lookup_key_value_from_file(path, pszProj, name,
sizeof(name)) > 0)
G_set_key_value("name", name, *projinfo);
else
G_set_key_value("name", pszProj, *projinfo);
}
/* the original hSRS1 is left as is, ok? */
hSRS = hSRS2;
use_proj_extension = 1;
}
G_free(proj4ext);
}
}
/* -------------------------------------------------------------------- */
/* Set cellhd for well known coordinate systems. */
/* -------------------------------------------------------------------- */
if (!OSRIsGeographic(hSRS) && !OSRIsProjected(hSRS))
goto default_to_xy;
if (cellhd) {
int bNorth;
if (OSRIsGeographic(hSRS)) {
cellhd->proj = PROJECTION_LL;
cellhd->zone = 0;
}
else if (OSRGetUTMZone(hSRS, &bNorth) != 0) {
cellhd->proj = PROJECTION_UTM;
cellhd->zone = OSRGetUTMZone(hSRS, &bNorth);
if (!bNorth)
cellhd->zone *= -1;
}
else {
cellhd->proj = PROJECTION_OTHER;
cellhd->zone = 0;
}
}
/* -------------------------------------------------------------------- */
/* Get the coordinate system definition in PROJ.4 format. */
/* -------------------------------------------------------------------- */
if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE)
goto default_to_xy;
/* -------------------------------------------------------------------- */
/* Parse the PROJ.4 string into key/value pairs. Do a bit of */
/* extra work to "GRASSify" the result. */
/* -------------------------------------------------------------------- */
temp_projinfo = G_create_key_value();
temp_projinfo_ext = G_create_key_value();
/* Create "local" copy of proj4 string so we can modify and free it
* using GRASS functions */
pszRemaining = G_store(pszProj4);
CPLFree(pszProj4);
pszProj4 = pszRemaining;
while ((pszRemaining = strstr(pszRemaining, "+")) != NULL) {
char *pszToken, *pszValue;
pszRemaining++;
/* Advance pszRemaining to end of this token[=value] pair */
pszToken = pszRemaining;
while (*pszRemaining != ' ' && *pszRemaining != '\0')
pszRemaining++;
if (*pszRemaining == ' ') {
*pszRemaining = '\0';
pszRemaining++;
}
/* parse token, and value */
if (strstr(pszToken, "=") != NULL) {
pszValue = strstr(pszToken, "=");
*pszValue = '\0';
pszValue++;
}
else
pszValue = "defined";
/* projection name */
if (G_strcasecmp(pszToken, "proj") == 0) {
/* The ll projection is known as longlat in PROJ.4 */
if (G_strcasecmp(pszValue, "longlat") == 0)
pszValue = "ll";
pszProj = pszValue;
}
/* Ellipsoid and datum handled separately below */
if (G_strcasecmp(pszToken, "ellps") == 0 ||
G_strcasecmp(pszToken, "a") == 0 ||
G_strcasecmp(pszToken, "b") == 0 ||
G_strcasecmp(pszToken, "es") == 0 ||
G_strcasecmp(pszToken, "rf") == 0 ||
G_strcasecmp(pszToken, "datum") == 0) {
G_set_key_value(pszToken, pszValue, temp_projinfo_ext);
continue;
}
/* We will handle units separately */
if (G_strcasecmp(pszToken, "to_meter") == 0)
continue;
if (G_strcasecmp(pszToken, "units") == 0) {
proj4_unit = G_store(pszValue);
continue;
}
G_set_key_value(pszToken, pszValue, temp_projinfo);
}
if (!pszProj)
G_warning(_("No projection name! Projection parameters likely to be "
"meaningless."));
/* -------------------------------------------------------------------- */
/* Derive the user name for the coordinate system. */
/* -------------------------------------------------------------------- */
if (!G_find_key_value("name", *projinfo)) {
pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
if (!pszProjCS)
pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
if (pszProjCS) {
G_set_key_value("name", pszProjCS, *projinfo);
}
else if (pszProj) {
char path[4095];
char name[80];
/* use name of the projection as name for the coordinate system */
snprintf(path, sizeof(path), "%s/etc/proj/projections",
G_gisbase());
if (G_lookup_key_value_from_file(path, pszProj, name,
sizeof(name)) > 0)
G_set_key_value("name", name, *projinfo);
else
G_set_key_value("name", pszProj, *projinfo);
}
}
/* -------------------------------------------------------------------- */
/* Find the GRASS datum name and choose parameters either */
/* interactively or not. */
/* -------------------------------------------------------------------- */
{
const char *pszDatumNameConst;
struct datum_list *list, *listhead;
char *dum1, *dum2, *pszDatumName;
int paramspresent = GPJ__get_datum_params(temp_projinfo, &dum1, &dum2);
if (!use_proj_extension)
pszDatumNameConst = OSRGetAttrValue(hSRS, "DATUM", 0);
else
pszDatumNameConst = G_find_key_value("datum", temp_projinfo_ext);
if (pszDatumNameConst) {
/* Need to make a new copy of the string so we don't mess
* around with the memory inside the OGRSpatialReferenceH? */
pszDatumName = G_store(pszDatumNameConst);
DatumNameMassage(&pszDatumName);
G_debug(3, "GPJ_osr_to_grass: pszDatumNameConst: <%s>",
pszDatumName);
list = listhead = read_datum_table();
while (list != NULL) {
if (G_strcasecmp(pszDatumName, list->longname) == 0) {
datum = G_store(list->name);
break;
}
list = list->next;
}
free_datum_list(listhead);
if (datum == NULL) {
if (paramspresent < 2)
/* Only give warning if no parameters present */
G_debug(1,
"Datum <%s> not recognised by GRASS and no "
"parameters found",
pszDatumName);
}
else {
G_set_key_value("datum", datum, *projinfo);
if (paramspresent < 2) {
/* If no datum parameters were imported from the OSR
* object then we should use the set specified by datumtrans
*/
char *params, *chosenparams = NULL;
int paramsets;
paramsets =
GPJ_get_default_datum_params_by_name(datum, &params);
if (paramsets < 0)
G_debug(1,
"Datum <%s> apparently recognised by GRASS but "
"no parameters found. "
"You may want to look into this.",
datum);
else if (datumtrans > paramsets) {
G_debug(
1,
"Invalid transformation number %d; valid range is "
"1 to %d. "
"Leaving datum transform parameters unspecified.",
datumtrans, paramsets);
datumtrans = 0;
}
if (paramsets > 0) {
struct gpj_datum_transform_list *tlist, *old;
tlist = GPJ_get_datum_transform_by_name(datum);
if (tlist != NULL) {
do {
if (tlist->count == datumtrans)
chosenparams = G_store(tlist->params);
old = tlist;
tlist = tlist->next;
GPJ_free_datum_transform(old);
} while (tlist != NULL);
}
}
if (chosenparams != NULL) {
char *paramkey, *paramvalue;
paramkey = strtok(chosenparams, "=");
paramvalue = chosenparams + strlen(paramkey) + 1;
G_set_key_value(paramkey, paramvalue, *projinfo);
G_free(chosenparams);
}
if (paramsets > 0)
G_free(params);
}
}
G_free(pszDatumName);
}
}
/* -------------------------------------------------------------------- */
/* Determine an appropriate GRASS ellipsoid name if possible, or */
/* else just put a and es values into PROJ_INFO */
/* -------------------------------------------------------------------- */
if ((datum != NULL) && (GPJ_get_datum_by_name(datum, &dstruct) > 0)) {
/* Use ellps name associated with datum */
G_set_key_value("ellps", dstruct.ellps, *projinfo);
GPJ_free_datum(&dstruct);
G_free(datum);
}
else if (!use_proj_extension) {
/* If we can't determine the ellipsoid from the datum, derive it
* directly from "SPHEROID" parameters in WKT */
const char *pszSemiMajor = OSRGetAttrValue(hSRS, "SPHEROID", 1);
const char *pszInvFlat = OSRGetAttrValue(hSRS, "SPHEROID", 2);
if (pszSemiMajor != NULL && pszInvFlat != NULL) {
char *ellps = NULL;
struct ellps_list *list, *listhead;
double a = atof(pszSemiMajor), invflat = atof(pszInvFlat), flat;
double es;
/* Allow for incorrect WKT describing a sphere where InvFlat
* is given as 0 rather than inf */
if (invflat > 0)
flat = 1 / invflat;
else
flat = 0;
es = flat * (2.0 - flat);
list = listhead = read_ellipsoid_table(0);
while (list != NULL) {
/* Try and match a and es against GRASS defined ellipsoids;
* accept first one that matches. These numbers were found
* by trial and error and could be fine-tuned, or possibly
* a direct comparison of IEEE floating point values used. */
if ((a == list->a || fabs(a - list->a) < 0.1 ||
fabs(1 - a / list->a) < 0.0000001) &&
((es == 0 && list->es == 0) ||
/* Special case for sphere */
(invflat == list->rf ||
fabs(invflat - list->rf) < 0.0000001))) {
ellps = G_store(list->name);
break;
}
list = list->next;
}
if (listhead != NULL)
free_ellps_list(listhead);
if (ellps == NULL) {
/* If we weren't able to find a matching ellps name, set
* a and es values directly from WKT-derived data */
char es_str[100];
G_set_key_value("a", (char *)pszSemiMajor, *projinfo);
snprintf(es_str, sizeof(es_str), "%.16g", es);
G_set_key_value("es", es_str, *projinfo);
}
else {
/* else specify the GRASS ellps name for readability */
G_set_key_value("ellps", ellps, *projinfo);
G_free(ellps);
}
}
}
else if (use_proj_extension) {
double a, es, rf;
if (GPJ__get_ellipsoid_params(temp_projinfo_ext, &a, &es, &rf)) {
char parmstr[100];
snprintf(parmstr, sizeof(parmstr), "%.16g", a);
G_set_key_value("a", parmstr, *projinfo);
snprintf(parmstr, sizeof(parmstr), "%.16g", es);
G_set_key_value("es", parmstr, *projinfo);
}
}
/* -------------------------------------------------------------------- */
/* Finally append the detailed projection parameters to the end */
/* -------------------------------------------------------------------- */
{
int i;
for (i = 0; i < temp_projinfo->nitems; i++)
G_set_key_value(temp_projinfo->key[i], temp_projinfo->value[i],
*projinfo);
G_free_key_value(temp_projinfo);
}
G_free_key_value(temp_projinfo_ext);
G_free(pszProj4);
/* -------------------------------------------------------------------- */
/* Set the linear units. */
/* -------------------------------------------------------------------- */
*projunits = G_create_key_value();
if (OSRIsGeographic(hSRS)) {
/* We assume degrees ... someday we will be wrong! */
G_set_key_value("unit", "degree", *projunits);
G_set_key_value("units", "degrees", *projunits);
G_set_key_value("meters", "1.0", *projunits);
}
else {
char szFormatBuf[256];
char *pszUnitsName = NULL;
double dfToMeters;
char *pszUnitsPlural, *pszStringEnd;
dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
/* the unit name can be arbitrary: the following can be the same
* us-ft (proj.4 keyword)
* U.S. Surveyor's Foot (proj.4 name)
* US survey foot (WKT)
* Foot_US (WKT)
*/
/* Workaround for the most obvious case when unit name is unknown */
if ((G_strcasecmp(pszUnitsName, "unknown") == 0) && (dfToMeters == 1.))
G_asprintf(&pszUnitsName, "meter");
if ((G_strcasecmp(pszUnitsName, "metre") == 0))
G_asprintf(&pszUnitsName, "meter");
if ((G_strcasecmp(pszUnitsName, "kilometre") == 0))
G_asprintf(&pszUnitsName, "kilometer");
if (dfToMeters != 1. && proj4_unit) {
int i;
i = 0;
while (gpj_units[i].id != NULL) {
if (strcmp(proj4_unit, gpj_units[i].id) == 0) {
G_asprintf(&pszUnitsName, "%s", gpj_units[i].name);
break;
}
i++;
}
}
G_set_key_value("unit", pszUnitsName, *projunits);
/* Attempt at plural formation (WKT format doesn't store plural
* form of unit name) */
pszUnitsPlural = G_malloc(strlen(pszUnitsName) + 3);
strcpy(pszUnitsPlural, pszUnitsName);
pszStringEnd = pszUnitsPlural + strlen(pszUnitsPlural) - 4;
if (G_strcasecmp(pszStringEnd, "foot") == 0) {
/* Special case for foot - change two o's to e's */
pszStringEnd[1] = 'e';
pszStringEnd[2] = 'e';
}
else if (G_strcasecmp(pszStringEnd, "inch") == 0) {
/* Special case for inch - add es */
pszStringEnd[4] = 'e';
pszStringEnd[5] = 's';
pszStringEnd[6] = '\0';
}
else {
/* For anything else add an s at the end */
pszStringEnd[4] = 's';
pszStringEnd[5] = '\0';
}
G_set_key_value("units", pszUnitsPlural, *projunits);
G_free(pszUnitsPlural);
snprintf(szFormatBuf, sizeof(szFormatBuf), "%.16g", dfToMeters);
G_set_key_value("meters", szFormatBuf, *projunits);
}
if (hSRS != hSRS1)
OSRDestroySpatialReference(hSRS);
return 2;
/* -------------------------------------------------------------------- */
/* Fallback to returning an ungeoreferenced definition. */
/* -------------------------------------------------------------------- */
default_to_xy:
if (cellhd != NULL) {
cellhd->proj = PROJECTION_XY;
cellhd->zone = 0;
}
if (*projinfo)
G_free_key_value(*projinfo);
*projinfo = NULL;
*projunits = NULL;
if (hSRS != NULL && hSRS != hSRS1)
OSRDestroySpatialReference(hSRS);
return 1;
}
/*!
* \brief Converts a WKT projection description to a GRASS co-ordinate system.
*
* \param cellhd Pointer to a GRASS Cell_head structure that will have its
* projection-related members populated with appropriate
* values \param projinfo Pointer to a pointer which will have a GRASS
* Key_Value structure allocated containing a set of GRASS PROJ_INFO values
* \param projunits Pointer to a pointer which will have a GRASS Key_Value
* structure allocated containing a set of GRASS PROJ_UNITS
* values \param wkt Well-known Text (WKT) description of the
* co-ordinate system to be converted \param datumtrans Index number of datum
* parameter set to use, 0 to leave unspecified
*
* \return 2 if a projected or lat/long co-ordinate system has been
* defined
* \return 1 if an unreferenced XY co-ordinate system has
* been defined
* \return -1 on error
*/
int GPJ_wkt_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
struct Key_Value **projunits, const char *wkt,
int datumtrans)
{
int retval;
if (wkt == NULL)
retval =
GPJ_osr_to_grass(cellhd, projinfo, projunits, NULL, datumtrans);
else {
OGRSpatialReferenceH hSRS;
/* Set finder function for locating OGR csv co-ordinate system tables */
/* SetCSVFilenameHook(GPJ_set_csv_loc); */
hSRS = OSRNewSpatialReference(wkt);
retval =
GPJ_osr_to_grass(cellhd, projinfo, projunits, hSRS, datumtrans);
OSRDestroySpatialReference(hSRS);
}
return retval;
}
/* GPJ_set_csv_loc()
* 'finder function' for use with OGR SetCSVFilenameHook() function */
const char *GPJ_set_csv_loc(const char *name)
{
const char *gisbase = G_gisbase();
static char *buf = NULL;
if (buf != NULL)
G_free(buf);
G_asprintf(&buf, "%s%s/%s", gisbase, CSVDIR, name);
return buf;
}
/* The list below is only for files that use a non-standard name for a
* datum that is already supported in GRASS. The number of entries must be even;
* they are all in pairs. The first one in the pair is the non-standard name;
* the second is the GRASS/GDAL name. If a name appears more than once (as for
* European_Terrestrial_Reference_System_1989) then it means there was more
* than one non-standard name for it that needs to be accounted for.
*
* N.B. The order of these pairs is different from that in
* ogr/ogrfromepsg.cpp in the GDAL source tree! GRASS uses the EPSG
* names in its WKT representation except WGS_1984 and WGS_1972 as
* these shortened versions seem to be standard.
* Below order:
* the equivalent name comes first in the pair, and
* the EPSG name (as used in the GRASS datum.table file) comes second.
*
* The datum parameters are stored in
* ../gis/datum.table # 3 parameters
* ../gis/datumtransform.table # 7 parameters (requires entry in datum.table)
*
* Hint: use GDAL's "testepsg" to identify the canonical name, e.g.
* testepsg epsg:4674
*/
static const char *papszDatumEquiv[] = {
"Militar_Geographische_Institute",
"Militar_Geographische_Institut",
"World_Geodetic_System_1984",
"WGS_1984",
"World_Geodetic_System_1972",
"WGS_1972",
"European_Terrestrial_Reference_System_89",
"European_Terrestrial_Reference_System_1989",
"European_Reference_System_1989",
"European_Terrestrial_Reference_System_1989",
"ETRS_1989",
"European_Terrestrial_Reference_System_1989",
"ETRS89",
"European_Terrestrial_Reference_System_1989",
"ETRF_1989",
"European_Terrestrial_Reference_System_1989",
"NZGD_2000",
"New_Zealand_Geodetic_Datum_2000",
"Monte_Mario_Rome",
"Monte_Mario",
"MONTROME",
"Monte_Mario",
"Campo_Inchauspe_1969",
"Campo_Inchauspe",
"S_JTSK",
"System_Jednotne_Trigonometricke_Site_Katastralni",
"S_JTSK_Ferro",
"Militar_Geographische_Institut",
"Potsdam_Datum_83",
"Deutsches_Hauptdreiecksnetz",
"Rauenberg_Datum_83",
"Deutsches_Hauptdreiecksnetz",
"South_American_1969",
"South_American_Datum_1969",
"International_Terrestrial_Reference_Frame_1992",
"ITRF92",
"ITRF_1992",
"ITRF92",
NULL};
/************************************************************************/
/* OGREPSGDatumNameMassage() */
/* */
/* Massage an EPSG datum name into WMT format. Also transform */
/* specific exception cases into WKT versions. */
/************************************************************************/
static void DatumNameMassage(char **ppszDatum)
{
int i, j;
char *pszDatum = *ppszDatum;
G_debug(3, "DatumNameMassage: Raw string found <%s>", (char *)pszDatum);
/* -------------------------------------------------------------------- */
/* Translate non-alphanumeric values to underscores. */
/* -------------------------------------------------------------------- */
for (i = 0; pszDatum[i] != '\0'; i++) {
if (!(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z') &&
!(pszDatum[i] >= 'a' && pszDatum[i] <= 'z') &&
!(pszDatum[i] >= '0' && pszDatum[i] <= '9')) {
pszDatum[i] = '_';
}
}
/* -------------------------------------------------------------------- */
/* Remove repeated and trailing underscores. */
/* -------------------------------------------------------------------- */
for (i = 1, j = 0; pszDatum[i] != '\0'; i++) {
if (pszDatum[j] == '_' && pszDatum[i] == '_')
continue;
pszDatum[++j] = pszDatum[i];
}
if (pszDatum[j] == '_')
pszDatum[j] = '\0';
else
pszDatum[j + 1] = '\0';
/* -------------------------------------------------------------------- */
/* Search for datum equivalences. Specific massaged names get */
/* mapped to OpenGIS specified names. */
/* -------------------------------------------------------------------- */
G_debug(3, "DatumNameMassage: Search for datum equivalences of <%s>",
(char *)pszDatum);
for (i = 0; papszDatumEquiv[i] != NULL; i += 2) {
if (EQUAL(*ppszDatum, papszDatumEquiv[i])) {
G_free(*ppszDatum);
*ppszDatum = G_store(papszDatumEquiv[i + 1]);
break;
}
}
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
G_asprintf
int G_asprintf(char **out, const char *fmt,...)
Definition asprintf.c:68
NULL
#define NULL
Definition ccmath.h:32
CSVDIR
#define CSVDIR
Definition convert.c:29
GPJ_grass_to_wkt
char * GPJ_grass_to_wkt(const struct Key_Value *proj_info, const struct Key_Value *proj_units, int esri_style, int prettify)
Converts a GRASS co-ordinate system representation to WKT style.
Definition convert.c:112
GPJ_grass_to_osr2
OGRSpatialReferenceH GPJ_grass_to_osr2(const struct Key_Value *proj_info, const struct Key_Value *proj_units, const struct Key_Value *proj_epsg)
Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object. EPSG code is preferred if avai...
Definition convert.c:347
GPJ_osr_to_grass
int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, OGRSpatialReferenceH hSRS1, int datumtrans)
Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
Definition convert.c:413
GPJ_wkt_to_grass
int GPJ_wkt_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, const char *wkt, int datumtrans)
Converts a WKT projection description to a GRASS co-ordinate system.
Definition convert.c:971
GPJ_grass_to_wkt2
char * GPJ_grass_to_wkt2(const struct Key_Value *proj_info, const struct Key_Value *proj_units, const struct Key_Value *proj_epsg, int esri_style, int prettify)
Converts a GRASS co-ordinate system representation to WKT style. EPSG code is preferred if available.
Definition convert.c:143
GPJ_grass_to_osr
OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value *proj_info, const struct Key_Value *proj_units)
Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
Definition convert.c:160
GPJ_set_csv_loc
const char * GPJ_set_csv_loc(const char *name)
Definition convert.c:998
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
GPJ_get_ellipsoid_by_name
int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
Looks up ellipsoid in ellipsoid table and returns the a, e2 parameters for the ellipsoid.
Definition ellipse.c:159
read_ellipsoid_table
struct ellps_list * read_ellipsoid_table(int fatal)
Definition ellipse.c:223
GPJ__get_ellipsoid_params
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a, double *e2, double *rf)
Get the ellipsoid parameters from proj keys structure.
Definition ellipse.c:74
free_ellps_list
void free_ellps_list(struct ellps_list *elist)
Definition ellipse.c:309
GPJ_free_ellps
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
Definition ellipse.c:302
pj_get_kv
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys, const struct Key_Value *in_units_keys)
Create a pj_info struct Co-ordinate System definition from a set of PROJ_INFO / PROJ_UNITS-style key-...
Definition get_proj.c:60
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
G_gisbase
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition gisbase.c:39
EQUAL
#define EQUAL(a, b)
Definition gsdrape.c:65
G_free_key_value
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition key_value1.c:104
G_set_key_value
void G_set_key_value(const char *key, const char *value, struct Key_Value *kv)
Set value for given key.
Definition key_value1.c:39
G_find_key_value
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive).
Definition key_value1.c:85
G_create_key_value
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition key_value1.c:23
G_lookup_key_value_from_file
int G_lookup_key_value_from_file(const char *file, const char *key, char value[], int n)
Look up for key in file.
Definition key_value4.c:48
name
const char * name
Definition named_colr.c:6
strcpy
#define strcpy
Definition parson.c:66
GPJ_get_default_datum_params_by_name
int GPJ_get_default_datum_params_by_name(const char *name, char **params)
"Last resort" function to retrieve a "default" set of datum parameters for a datum (N....
Definition proj/datum.c:86
GPJ_free_datum_transform
void GPJ_free_datum_transform(struct gpj_datum_transform_list *item)
Free the memory used by a gpj_datum_transform_list struct.
Definition proj/datum.c:323
GPJ_free_datum
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
Definition proj/datum.c:396
GPJ_get_datum_transform_by_name
struct gpj_datum_transform_list * GPJ_get_datum_transform_by_name(const char *inputname)
Internal function to find all possible sets of transformation parameters for a particular datum.
Definition proj/datum.c:237
free_datum_list
void free_datum_list(struct datum_list *dstruct)
Free the memory used by a datum_list linked list structure.
Definition proj/datum.c:410
GPJ__get_datum_params
int GPJ__get_datum_params(const struct Key_Value *projinfo, char **datumname, char **params)
Extract the datum transformation-related parameters from a set of general PROJ_INFO parameters.
Definition proj/datum.c:173
GPJ_get_datum_by_name
int GPJ_get_datum_by_name(const char *name, struct gpj_datum *dstruct)
Look up a string in datum.table file to see if it is a valid datum name and if so place its informati...
Definition proj/datum.c:38
read_datum_table
struct datum_list * read_datum_table(void)
Read the current GRASS datum.table from disk and store in memory.
Definition proj/datum.c:342
list
struct list * list
Definition read_list.c:24
G_strcasecmp
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower).
Definition strings.c:47
G_store
char * G_store(const char *s)
Copy string to allocated memory.
Definition strings.c:87
path
Definition path.h:15
G_free_tokens
void G_free_tokens(char **tokens)
Free memory allocated to tokens.
Definition token.c:197
G_tokenize
char ** G_tokenize(const char *buf, const char *delim)
Tokenize string.
Definition token.c:47
G_number_of_tokens
int G_number_of_tokens(char **tokens)
Return number of tokens.
Definition token.c:178