16 #include <grass/config.h> 22 #include <grass/gis.h> 23 #include <grass/gprojects.h> 24 #include <grass/glocale.h> 28 #include "local_proto.h" 31 #define CSVDIR "/etc/proj/ogr_csv" 33 static void DatumNameMassage(
char **);
46 {
"km",
"1000.",
"Kilometer", 1000.0},
47 {
"m",
"1.",
"Meter", 1.0},
48 {
"dm",
"1/10",
"Decimeter", 0.1},
49 {
"cm",
"1/100",
"Centimeter", 0.01},
50 {
"mm",
"1/1000",
"Millimeter", 0.001},
51 {
"kmi",
"1852.0",
"International Nautical Mile", 1852.0},
52 {
"in",
"0.0254",
"International Inch", 0.0254},
53 {
"ft",
"0.3048",
"International Foot", 0.3048},
54 {
"yd",
"0.9144",
"International Yard", 0.9144},
55 {
"mi",
"1609.344",
"International Statute Mile", 1609.344},
56 {
"fath",
"1.8288",
"International Fathom", 1.8288},
57 {
"ch",
"20.1168",
"International Chain", 20.1168},
58 {
"link",
"0.201168",
"International Link", 0.201168},
59 {
"us-in",
"1./39.37",
"U.S. Surveyor's Inch", 0.0254},
60 {
"us-ft",
"0.304800609601219",
"U.S. Surveyor's Foot", 0.304800609601219},
61 {
"us-yd",
"0.914401828803658",
"U.S. Surveyor's Yard", 0.914401828803658},
62 {
"us-ch",
"20.11684023368047",
"U.S. Surveyor's Chain", 20.11684023368047},
63 {
"us-mi",
"1609.347218694437",
"U.S. Surveyor's Statute Mile", 1609.347218694437},
64 {
"ind-yd",
"0.91439523",
"Indian Yard", 0.91439523},
65 {
"ind-ft",
"0.30479841",
"Indian Foot", 0.30479841},
66 {
"ind-ch",
"20.11669506",
"Indian Chain", 20.11669506},
70 static char *grass_to_wkt(
const struct Key_Value *proj_info,
71 const struct Key_Value *proj_units,
72 const struct Key_Value *proj_epsg,
73 int esri_style,
int prettify)
76 OGRSpatialReferenceH hSRS;
77 char *wkt, *local_wkt;
88 OSRExportToPrettyWkt(hSRS, &wkt, 0);
90 OSRExportToWkt(hSRS, &wkt);
94 OSRDestroySpatialReference(hSRS);
98 G_warning(_(
"GRASS is not compiled with OGR support"));
122 const struct Key_Value *proj_units,
123 int esri_style,
int prettify)
125 return grass_to_wkt(proj_info, proj_units,
NULL, esri_style, prettify);
153 const struct Key_Value *proj_units,
154 const struct Key_Value *proj_epsg,
155 int esri_style,
int prettify)
157 return grass_to_wkt(proj_info, proj_units, proj_epsg, esri_style, prettify);
171 const struct Key_Value * proj_units)
173 struct pj_info pjinfo;
174 char *proj4, *proj4mod, *wkt, *modwkt, *startmod, *lastpart;
175 OGRSpatialReferenceH hSRS, hSRS2;
177 struct gpj_datum dstruct;
178 struct gpj_ellps estruct;
180 const char *ellpskv, *unit, *unfact;
181 char *ellps, *ellpslong, *datum, *params, *towgs84, *datumlongname,
183 const char *sysname, *osrunit, *osrunfact;
187 if ((proj_info ==
NULL) || (proj_units ==
NULL))
190 hSRS = OSRNewSpatialReference(
NULL);
193 if (
pj_get_kv(&pjinfo, proj_info, proj_units) < 0) {
194 G_warning(_(
"Unable parse GRASS PROJ_INFO file"));
200 if ((proj4 = pjinfo.def) ==
NULL) {
201 G_warning(_(
"Unable get PROJ.4-style parameter string"));
205 proj_destroy(pjinfo.pj);
212 if (unfact !=
NULL && (strcmp(pjinfo.proj,
"ll") != 0))
213 G_asprintf(&proj4mod,
"%s +to_meter=%s", proj4, unfact);
218 if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
219 G_warning(_(
"OGR can't parse PROJ.4-style parameter string: " 220 "%s (OGR Error code was %d)"), proj4mod, errcode);
231 if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
232 G_warning(_(
"OGR can't get WKT-style parameter string " 233 "(OGR Error code was %d)"), errcode);
247 datumlongname =
G_store(
"unknown");
252 datumlongname =
G_store(dstruct.longname);
254 ellps =
G_store(dstruct.ellps);
257 G_debug(3,
"GPJ_grass_to_osr: datum: <%s>", datum);
260 ellpslong =
G_store(estruct.longname);
261 DatumNameMassage(&ellpslong);
267 startmod = strstr(wkt,
"GEOGCS");
268 lastpart = strstr(wkt,
"PRIMEM");
269 len = strlen(wkt) - strlen(startmod);
271 if (haveparams == 2) {
274 char *paramkey, *paramvalue;
276 paramkey = strtok(params,
"=");
277 paramvalue = params + strlen(paramkey) + 1;
279 G_asprintf(&towgs84,
",TOWGS84[%s]", paramvalue);
287 sysname = OSRGetAttrValue(hSRS,
"PROJCS", 0);
288 if (sysname ==
NULL) {
294 if ((strcmp(sysname,
"unnamed") == 0) &&
301 osrunit = OSRGetAttrValue(hSRS,
"UNIT", 0);
302 osrunfact = OSRGetAttrValue(hSRS,
"UNIT", 1);
308 double unfactf = atof(unfact);
312 startmod = strstr(lastpart, buff);
313 len = strlen(lastpart) - strlen(startmod);
314 lastpart[len] =
'\0';
319 G_asprintf(&end,
",UNIT[\"%s\",%.16g]]", unit, unfactf);
323 OSRDestroySpatialReference(hSRS);
325 "%sGEOGCS[\"%s\",DATUM[\"%s\",SPHEROID[\"%s\",%.16g,%.16g]%s],%s%s",
326 start, ellps, datumlongname, ellpslong, a, rf, towgs84,
328 hSRS2 = OSRNewSpatialReference(modwkt);
362 const struct Key_Value * proj_units,
363 const struct Key_Value * proj_epsg)
370 epsgcode = atoi(epsgstr);
375 OGRSpatialReferenceH hSRS;
377 hSRS = OSRNewSpatialReference(
NULL);
379 OSRImportFromEPSG(hSRS, epsgcode);
386 double df[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
391 df[i] = atof(tokens[i]);
394 OSRSetTOWGS84(hSRS, df[0], df[1], df[2], df[3], df[4], df[5], df[6]);
423 struct Key_Value **projunits, OGRSpatialReferenceH hSRS1,
426 struct Key_Value *temp_projinfo;
427 char *pszProj4 =
NULL, *pszRemaining;
428 char *pszProj =
NULL;
429 const char *pszProjCS =
NULL;
431 char *proj4_unit =
NULL;
432 struct gpj_datum dstruct;
434 OGRSpatialReferenceH hSRS;
449 OSRMorphFromESRI(hSRS);
454 ograttr = OSRGetAttrValue(hSRS,
"EXTENSION", 0);
455 if (ograttr && *ograttr && strcmp(ograttr,
"PROJ4") == 0) {
456 ograttr = OSRGetAttrValue(hSRS,
"EXTENSION", 1);
457 G_debug(3,
"proj4 extension:");
460 if (ograttr && *ograttr) {
462 OGRSpatialReferenceH hSRS2;
464 hSRS2 = OSRNewSpatialReference(
NULL);
468 if (OSRImportFromProj4(hSRS2, proj4ext) != OGRERR_NONE) {
469 G_warning(_(
"Updating spatial reference with embedded proj4 definition failed. " 470 "Proj4 definition: <%s>"), proj4ext);
471 OSRDestroySpatialReference(hSRS2);
476 G_warning(_(
"Updating spatial reference with embedded proj4 definition"));
481 pszProjCS = OSRGetAttrValue(hSRS,
"PROJCS", 0);
483 pszProjCS = OSRGetAttrValue(hSRS,
"GEOGCS", 0);
494 sprintf(path,
"%s/etc/proj/projections",
G_gisbase());
512 if (!OSRIsGeographic(hSRS) && !OSRIsProjected(hSRS))
518 if (OSRIsGeographic(hSRS)) {
519 cellhd->proj = PROJECTION_LL;
522 else if (OSRGetUTMZone(hSRS, &bNorth) != 0) {
523 cellhd->proj = PROJECTION_UTM;
524 cellhd->zone = OSRGetUTMZone(hSRS, &bNorth);
529 cellhd->proj = PROJECTION_OTHER;
537 if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE)
548 pszRemaining =
G_store(pszProj4);
550 pszProj4 = pszRemaining;
551 while ((pszRemaining = strstr(pszRemaining,
"+")) !=
NULL) {
552 char *pszToken, *pszValue;
557 pszToken = pszRemaining;
558 while (*pszRemaining !=
' ' && *pszRemaining !=
'\0')
561 if (*pszRemaining ==
' ') {
562 *pszRemaining =
'\0';
567 if (strstr(pszToken,
"=") !=
NULL) {
568 pszValue = strstr(pszToken,
"=");
573 pszValue =
"defined";
603 proj4_unit =
G_store(pszValue);
610 G_warning(_(
"No projection name! Projection parameters likely to be meaningless."));
616 pszProjCS = OSRGetAttrValue(hSRS,
"PROJCS", 0);
618 pszProjCS = OSRGetAttrValue(hSRS,
"GEOGCS", 0);
629 sprintf(path,
"%s/etc/proj/projections",
G_gisbase());
644 const char *pszDatumNameConst = OSRGetAttrValue(hSRS,
"DATUM", 0);
645 struct datum_list *
list, *listhead;
646 char *dum1, *dum2, *pszDatumName;
650 if (pszDatumNameConst) {
654 pszDatumName =
G_store(pszDatumNameConst);
655 DatumNameMassage(&pszDatumName);
656 G_debug(3,
"GPJ_osr_to_grass: pszDatumNameConst: <%s>", pszDatumName);
660 while (list !=
NULL) {
670 if (paramspresent < 2)
672 G_warning(_(
"Datum <%s> not recognised by GRASS and no parameters found"),
678 if (paramspresent < 2) {
681 char *params, *chosenparams =
NULL;
688 G_warning(_(
"Datum <%s> apparently recognised by GRASS but no parameters found. " 689 "You may want to look into this."), datum);
690 else if (datumtrans > paramsets) {
692 G_warning(_(
"Invalid transformation number %d; valid range is 1 to %d. " 693 "Leaving datum transform parameters unspecified."),
694 datumtrans, paramsets);
699 struct gpj_datum_transform_list *tlist, *old;
705 if (tlist->count == datumtrans)
706 chosenparams =
G_store(tlist->params);
710 }
while (tlist !=
NULL);
714 if (chosenparams !=
NULL) {
715 char *paramkey, *paramvalue;
717 paramkey = strtok(chosenparams,
"=");
718 paramvalue = chosenparams + strlen(paramkey) + 1;
746 const char *pszSemiMajor = OSRGetAttrValue(hSRS,
"SPHEROID", 1);
747 const char *pszInvFlat = OSRGetAttrValue(hSRS,
"SPHEROID", 2);
749 if (pszSemiMajor !=
NULL && pszInvFlat !=
NULL) {
751 struct ellps_list *
list, *listhead;
752 double a = atof(pszSemiMajor), invflat = atof(pszInvFlat), flat;
762 es = flat * (2.0 - flat);
766 while (list !=
NULL) {
771 if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) &&
772 ((es == 0 && list->es == 0) ||
774 (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
780 if (listhead !=
NULL)
790 sprintf(es_str,
"%.16g", es);
810 for (i = 0; i < temp_projinfo->nitems; i++)
812 temp_projinfo->value[i], *projinfo);
824 if (OSRIsGeographic(hSRS)) {
831 char szFormatBuf[256];
832 char *pszUnitsName =
NULL;
834 char *pszUnitsPlural, *pszStringEnd;
836 dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
855 if (dfToMeters != 1. && proj4_unit) {
859 while (gpj_units[i].
id !=
NULL) {
860 if (strcmp(proj4_unit, gpj_units[i].
id) == 0) {
872 pszUnitsPlural = G_malloc(strlen(pszUnitsName) + 3);
873 strcpy(pszUnitsPlural, pszUnitsName);
874 pszStringEnd = pszUnitsPlural + strlen(pszUnitsPlural) - 4;
877 pszStringEnd[1] =
'e';
878 pszStringEnd[2] =
'e';
882 pszStringEnd[4] =
'e';
883 pszStringEnd[5] =
's';
884 pszStringEnd[6] =
'\0';
888 pszStringEnd[4] =
's';
889 pszStringEnd[5] =
'\0';
895 sprintf(szFormatBuf,
"%.16g", dfToMeters);
901 OSRDestroySpatialReference(hSRS);
909 if (cellhd !=
NULL) {
910 cellhd->proj = PROJECTION_XY;
919 if (hSRS !=
NULL && hSRS != hSRS1)
920 OSRDestroySpatialReference(hSRS);
947 struct Key_Value **projunits,
const char *wkt,
957 OGRSpatialReferenceH hSRS;
962 hSRS = OSRNewSpatialReference(wkt);
965 OSRDestroySpatialReference(hSRS);
981 static char *buf =
NULL;
1015 static const char *papszDatumEquiv[] = {
1016 "Militar_Geographische_Institute",
1017 "Militar_Geographische_Institut",
1018 "World_Geodetic_System_1984",
1020 "World_Geodetic_System_1972",
1022 "European_Terrestrial_Reference_System_89",
1023 "European_Terrestrial_Reference_System_1989",
1024 "European_Reference_System_1989",
1025 "European_Terrestrial_Reference_System_1989",
1027 "European_Terrestrial_Reference_System_1989",
1029 "European_Terrestrial_Reference_System_1989",
1031 "European_Terrestrial_Reference_System_1989",
1033 "New_Zealand_Geodetic_Datum_2000",
1038 "Campo_Inchauspe_1969",
1041 "System_Jednotne_Trigonometricke_Site_Katastralni",
1043 "Militar_Geographische_Institut",
1045 "Deutsches_Hauptdreiecksnetz",
1046 "South_American_1969",
1047 "South_American_Datum_1969",
1061 static void DatumNameMassage(
char **ppszDatum)
1064 char *pszDatum = *ppszDatum;
1066 G_debug(3,
"DatumNameMassage: Raw string found <%s>", (
char *)pszDatum);
1070 for (i = 0; pszDatum[i] !=
'\0'; i++) {
1071 if (!(pszDatum[i] >=
'A' && pszDatum[i] <=
'Z')
1072 && !(pszDatum[i] >=
'a' && pszDatum[i] <=
'z')
1073 && !(pszDatum[i] >=
'0' && pszDatum[i] <=
'9')) {
1081 for (i = 1, j = 0; pszDatum[i] !=
'\0'; i++) {
1082 if (pszDatum[j] ==
'_' && pszDatum[i] ==
'_')
1085 pszDatum[++j] = pszDatum[i];
1087 if (pszDatum[j] ==
'_')
1090 pszDatum[j + 1] =
'\0';
1096 G_debug(3,
"DatumNameMassage: Search for datum equivalences of <%s>", (
char *)pszDatum);
1097 for (i = 0; papszDatumEquiv[i] !=
NULL; i += 2) {
1098 if (
EQUAL(*ppszDatum, papszDatumEquiv[i])) {
1100 *ppszDatum =
G_store(papszDatumEquiv[i + 1]);
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)
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
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...
void GPJ_free_datum_transform(struct gpj_datum_transform_list *item)
Free the memory used by a gpj_datum_transform_list struct.
struct ellps_list * read_ellipsoid_table(int fatal)
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.
char * G_store(const char *s)
Copy string to allocated memory.
int G_asprintf(char **out, const char *fmt,...)
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...
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...
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
struct gpj_units gpj_units[]
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.
int G_lookup_key_value_from_file(const char *file, const char *key, char value[], int n)
Look up for key in file.
int G_number_of_tokens(char **tokens)
Return number of tokens.
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)
int G_debug(int level, const char *msg,...)
Print debugging message.
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
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.
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...
void G_free_tokens(char **tokens)
Free memory allocated to tokens.
struct datum_list * read_datum_table(void)
Read the current GRASS datum.table from disk and store in memory.
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-...
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_tokenize(const char *buf, const char *delim)
Tokenize string.
void G_set_key_value(const char *key, const char *value, struct Key_Value *kv)
Set value for given key.
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.
void G_free(void *buf)
Free allocated memory.
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.
const char * G_gisbase(void)
Get full path name of the top level module directory.
void G_warning(const char *msg,...)
Print a warning message to stderr.
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
void free_ellps_list(struct ellps_list *elist)