16 #include <grass/config.h>
23 #include <grass/gis.h>
24 #include <grass/gprojects.h>
25 #include <grass/glocale.h>
30 #define CSVDIR "/etc/ogr_csv"
32 static void DatumNameMassage(
char **);
53 struct Key_Value *proj_units,
54 int esri_style,
int prettify)
56 OGRSpatialReferenceH hSRS;
57 char *wkt, *local_wkt;
68 OSRExportToPrettyWkt(hSRS, &wkt, 0);
70 OSRExportToWkt(hSRS, &wkt);
74 OSRDestroySpatialReference(hSRS);
89 struct Key_Value * proj_units)
91 struct pj_info pjinfo;
92 char *proj4, *proj4mod, *wkt, *modwkt, *startmod, *lastpart;
93 OGRSpatialReferenceH hSRS, hSRS2;
95 struct gpj_datum dstruct;
96 struct gpj_ellps estruct;
98 char *ellps, *ellpslong, *datum, *params, *towgs84, *datumlongname,
100 const char *sysname, *osrunit, *osrunfact;
104 if ((proj_info ==
NULL) || (proj_units ==
NULL))
107 hSRS = OSRNewSpatialReference(
NULL);
109 if (
pj_get_kv(&pjinfo, proj_info, proj_units) < 0) {
110 G_warning(_(
"Unable parse GRASS PROJ_INFO file"));
114 if ((proj4 = pj_get_def(pjinfo.pj, 0)) ==
NULL) {
115 G_warning(_(
"Unable get PROJ.4-style parameter string"));
121 if (unfact !=
NULL && (strcmp(pjinfo.proj,
"ll") != 0))
122 G_asprintf(&proj4mod,
"%s +to_meter=%s", proj4, unfact);
126 if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
127 G_warning(_(
"OGR can't parse PROJ.4-style parameter string: "
128 "%s (OGR Error code was %d)"), proj4mod, errcode);
132 if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
133 G_warning(_(
"OGR can't get WKT-style parameter string "
134 "(OGR Error code was %d)"), errcode);
148 datumlongname =
G_store(dstruct.longname);
150 ellps =
G_store(dstruct.ellps);
154 ellpslong =
G_store(estruct.longname);
155 DatumNameMassage(&ellpslong);
163 len = strlen(wkt) - strlen(startmod);
165 if (haveparams == 2) {
168 char *paramkey, *paramvalue;
170 paramkey = strtok(params,
"=");
171 paramvalue = params + strlen(paramkey) + 1;
173 G_asprintf(&towgs84,
",TOWGS84[%s]", paramvalue);
180 sysname = OSRGetAttrValue(hSRS,
"PROJCS", 0);
181 if (sysname ==
NULL) {
187 if ((strcmp(sysname,
"unnamed") == 0) &&
194 osrunit = OSRGetAttrValue(hSRS,
"UNIT", 0);
195 osrunfact = OSRGetAttrValue(hSRS,
"UNIT", 1);
200 double unfactf = atof(unfact);
204 startmod =
G_strstr(lastpart, buff);
205 len = strlen(lastpart) - strlen(startmod);
206 lastpart[len] =
'\0';
210 G_asprintf(&end,
",UNIT[\"%s\",%.16g]]", unit, unfactf);
216 "%sGEOGCS[\"%s\",DATUM[\"%s\",SPHEROID[\"%s\",%.16g,%.16g]%s],%s%s",
217 start, ellps, datumlongname, ellpslong, a, rf, towgs84,
220 hSRS2 = OSRNewSpatialReference(modwkt);
222 OSRDestroySpatialReference(hSRS);
226 if (proj4 != proj4mod)
259 struct Key_Value **projunits, OGRSpatialReferenceH hSRS,
262 struct Key_Value *temp_projinfo;
263 char *pszProj4 =
NULL, *pszRemaining;
264 char *pszProj =
NULL;
266 struct gpj_datum dstruct;
276 OSRMorphFromESRI(hSRS);
281 if (!OSRIsGeographic(hSRS) && !OSRIsProjected(hSRS))
287 if (OSRIsGeographic(hSRS)) {
288 cellhd->proj = PROJECTION_LL;
291 else if (OSRGetUTMZone(hSRS, &bNorth) != 0) {
292 cellhd->proj = PROJECTION_UTM;
293 cellhd->zone = OSRGetUTMZone(hSRS, &bNorth);
298 cellhd->proj = PROJECTION_OTHER;
306 if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE)
317 pszRemaining =
G_store(pszProj4);
319 pszProj4 = pszRemaining;
320 while ((pszRemaining = strstr(pszRemaining,
"+")) !=
NULL) {
321 char *pszToken, *pszValue;
326 pszToken = pszRemaining;
327 while (*pszRemaining !=
' ' && *pszRemaining !=
'\0')
330 if (*pszRemaining ==
' ') {
331 *pszRemaining =
'\0';
336 if (strstr(pszToken,
"=") !=
NULL) {
337 pszValue = strstr(pszToken,
"=");
342 pszValue =
"defined";
390 G_warning(_(
"No projection name! Projection parameters likely to be meaningless."));
399 const char *pszDatumNameConst = OSRGetAttrValue(hSRS,
"DATUM", 0);
401 char *dum1, *dum2, *pszDatumName;
405 if (pszDatumNameConst) {
409 pszDatumName =
G_store(pszDatumNameConst);
410 DatumNameMassage(&pszDatumName);
414 while (list !=
NULL) {
424 if (paramspresent < 2)
426 G_warning(_(
"Datum <%s> not recognised by GRASS and no parameters found"),
432 if (paramspresent < 2) {
435 char *params, *chosenparams =
NULL;
442 G_warning(_(
"Datum <%s> apparently recognised by GRASS but no parameters found. "
443 "You may want to look into this."), datum);
444 else if (datumtrans > paramsets) {
446 G_warning(_(
"Invalid transformation number %d; valid range is 1 to %d. "
447 "Leaving datum transform parameters unspecified."),
448 datumtrans, paramsets);
453 struct gpj_datum_transform_list *list, *
old;
459 if (list->count == datumtrans) {
460 chosenparams =
G_store(list->params);
466 }
while (list !=
NULL);
470 if (chosenparams !=
NULL) {
471 char *paramkey, *paramvalue;
473 paramkey = strtok(chosenparams,
"=");
474 paramvalue = chosenparams + strlen(paramkey) + 1;
501 const char *pszSemiMajor = OSRGetAttrValue(hSRS,
"SPHEROID", 1);
502 const char *pszInvFlat = OSRGetAttrValue(hSRS,
"SPHEROID", 2);
504 if (pszSemiMajor !=
NULL && pszInvFlat !=
NULL) {
507 double a = atof(pszSemiMajor), invflat = atof(pszInvFlat), flat;
517 es = flat * (2.0 - flat);
521 while (list !=
NULL) {
526 if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) && ((es == 0 && list->es == 0) ||
527 (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
533 if (listhead !=
NULL)
563 for (i = 0; i < temp_projinfo->nitems; i++)
565 temp_projinfo->value[i], *projinfo);
577 if (OSRIsGeographic(hSRS)) {
584 char szFormatBuf[256];
585 char *pszUnitsName =
NULL;
587 char *pszUnitsPlural, *pszStringEnd;
589 dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
600 pszUnitsPlural = G_malloc(strlen(pszUnitsName) + 3);
601 strcpy(pszUnitsPlural, pszUnitsName);
602 pszStringEnd = pszUnitsPlural + strlen(pszUnitsPlural) - 4;
605 pszStringEnd[1] =
'e';
606 pszStringEnd[2] =
'e';
610 pszStringEnd[4] =
'e';
611 pszStringEnd[5] =
's';
612 pszStringEnd[6] =
'\0';
616 pszStringEnd[4] =
's';
617 pszStringEnd[5] =
'\0';
623 sprintf(szFormatBuf,
"%.16g", dfToMeters);
634 if (cellhd !=
NULL) {
635 cellhd->proj = PROJECTION_XY;
666 struct Key_Value **projunits,
const char *wkt,
675 OGRSpatialReferenceH hSRS;
680 hSRS = OSRNewSpatialReference(wkt);
683 OSRDestroySpatialReference(hSRS);
720 static const char *papszDatumEquiv[] = {
721 "Militar_Geographische_Institute",
722 "Militar_Geographische_Institut",
723 "World_Geodetic_System_1984",
725 "World_Geodetic_System_1972",
727 "European_Terrestrial_Reference_System_89",
728 "European_Terrestrial_Reference_System_1989",
729 "European_Reference_System_1989",
730 "European_Terrestrial_Reference_System_1989",
732 "European_Terrestrial_Reference_System_1989",
734 "European_Terrestrial_Reference_System_1989",
736 "European_Terrestrial_Reference_System_1989",
738 "New_Zealand_Geodetic_Datum_2000",
743 "Campo_Inchauspe_1969",
746 "Militar_Geographische_Institut",
748 "Deutsches_Hauptdreiecksnetz",
749 "South_American_1969",
750 "South_American_Datum_1969",
764 static void DatumNameMassage(
char **ppszDatum)
767 char *pszDatum = *ppszDatum;
772 for (i = 0; pszDatum[i] !=
'\0'; i++) {
773 if (!(pszDatum[i] >=
'A' && pszDatum[i] <=
'Z')
774 && !(pszDatum[i] >=
'a' && pszDatum[i] <=
'z')
775 && !(pszDatum[i] >=
'0' && pszDatum[i] <=
'9')) {
783 for (i = 1, j = 0; pszDatum[i] !=
'\0'; i++) {
784 if (pszDatum[j] ==
'_' && pszDatum[i] ==
'_')
787 pszDatum[++j] = pszDatum[i];
789 if (pszDatum[j] ==
'_')
792 pszDatum[j + 1] =
'\0';
798 for (i = 0; papszDatumEquiv[i] !=
NULL; i += 2) {
799 if (
EQUAL(*ppszDatum, papszDatumEquiv[i])) {
801 *ppszDatum =
G_store(papszDatumEquiv[i + 1]);
char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key.
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.
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
sprintf(buf2,"%s", G3D_CATS_ELEMENT)
void G_free(void *buf)
Free allocated memory.
int pj_get_kv(struct pj_info *info, struct Key_Value *in_proj_keys, 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-...
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
char * G_store(const char *s)
Copy string to allocated memory.
struct ellps_list * read_ellipsoid_table(int fatal)
int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, OGRSpatialReferenceH hSRS, int datumtrans)
Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
int G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
int G_asprintf(char **out, const char *fmt,...)
char * GPJ_grass_to_wkt(struct Key_Value *proj_info, struct Key_Value *proj_units, int esri_style, int prettify)
Converts a GRASS co-ordinate system representation to WKT style.
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...
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...
int G_lookup_key_value_from_file(const char *file, const char *key, char value[], int n)
Look up for key in file.
OGRSpatialReferenceH GPJ_grass_to_osr(struct Key_Value *proj_info, struct Key_Value *proj_units)
Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
int GPJ__get_ellipsoid_params(struct Key_Value *proj_keys, double *a, double *e2, double *rf)
void free_datum_list(struct datum_list *dstruct)
Free the memory used by a datum_list linked list structure.
const char * GPJ_set_csv_loc(const char *name)
void GPJ_free_ellps(struct gpj_ellps *estruct)
struct datum_list * read_datum_table(void)
Read the current GRASS datum.table from disk and store in memory.
char buf[GNAME_MAX+sizeof(G3D_DIRECTORY)+2]
int GPJ__get_datum_params(struct Key_Value *projinfo, char **datumname, char **params)
Extract the datum transformation-related parameters from a set of general PROJ_INFO parameters...
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
char * G_strstr(const char *mainString, const char *subString)
Finds the first occurrence of the character C in the null-terminated string beginning at mainString...
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...
char * G_gisbase(void)
top level module directory
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 ...
int G_set_key_value(const char *key, const char *value, struct Key_Value *kv)
Set value for given key.
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
void free_ellps_list(struct ellps_list *elist)