1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /* Program: gmt_ogrproj.c
18 * Purpose: routine to do point coordinate transformations indirectly with the proj.4 lib
19 *
20 * Author: Joaquim Luis
21 * Date: 17-Aug-2017
22 */
23
24 #include <gdal.h>
25 #include <ogr_srs_api.h>
26
gmt_OGRCoordinateTransformation(struct GMT_CTRL * GMT,const char * pSrcSRS,const char * pDstSRS)27 OGRCoordinateTransformationH gmt_OGRCoordinateTransformation(struct GMT_CTRL *GMT, const char *pSrcSRS, const char *pDstSRS) {
28 /* pSrcSRS and pDstSRS are pointers to strings defining the Source and Destination Referencing
29 System. The SRS can be a +proj Proj.4 string, a WKT, a EPSG:n code or a filename with a WKT (?).
30
31 The caller to this function is responsible to free the GDAL object created here with a call to
32 OCTDestroyCoordinateTransformation(hCT);
33 */
34 OGRSpatialReferenceH hSrcSRS, hDstSRS;
35 OGRCoordinateTransformationH hCT;
36
37 /* ------------------ Set the Source projection ----------------------------- */
38 hSrcSRS = OSRNewSpatialReference(NULL);
39 if (OSRSetFromUserInput(hSrcSRS, pSrcSRS) != OGRERR_NONE) {
40 GMT_Report (GMT->parent, GMT_MSG_ERROR, "OGRPROJ: Translating source SRS failed.\n%s\n", pSrcSRS);
41 return NULL;
42 }
43 #if GDAL_VERSION_MAJOR >= 3
44 OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER); /* Set the data axis to CRS axis mapping strategy. */
45 #endif
46 /* ------------------- Set the Target projection ---------------------------- */
47 CPLErrorReset();
48 hDstSRS = OSRNewSpatialReference(NULL);
49 if (OSRSetFromUserInput(hDstSRS, pDstSRS) != OGRERR_NONE) {
50 GMT_Report (GMT->parent, GMT_MSG_ERROR, "OGRPROJ: Translating target SRS failed.\n%s\n", pDstSRS);
51 OSRDestroySpatialReference(hSrcSRS); /* It was just created above */
52 return NULL;
53 }
54 #if GDAL_VERSION_MAJOR >= 3
55 OSRSetAxisMappingStrategy(hDstSRS, OAMS_TRADITIONAL_GIS_ORDER); /* Set the data axis to CRS axis mapping strategy. */
56 #endif
57 /* -------------------------------------------------------------------------- */
58
59 hCT = OCTNewCoordinateTransformation(hSrcSRS, hDstSRS);
60 if (hCT == NULL) {
61 char *pszSrcWKT = NULL, *pszDstWKT = NULL;
62 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to create coordinate transformation between the following\n"
63 "coordinate systems. This may be because they are not transformable,\n"
64 "or because projection services (PROJ.4 DLL/.so) could not be loaded.\n");
65 OSRExportToPrettyWkt(hSrcSRS, &pszSrcWKT, FALSE);
66 OSRExportToPrettyWkt(hDstSRS, &pszDstWKT, FALSE);
67 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Source:\n\n%s\n\n%s\n\n", pszSrcWKT, pszDstWKT);
68 CPLFree(pszSrcWKT); CPLFree(pszDstWKT);
69 }
70 OSRDestroySpatialReference(hSrcSRS); OSRDestroySpatialReference(hDstSRS);
71 return hCT;
72 }
73
gmt_ogrproj(struct GMT_CTRL * GMT,char * pSrcSRS,char * pDstSRS,int n_pts,double * xi,double * yi,double * zi,bool insitu,double * xo,double * yo,double * zo)74 int gmt_ogrproj(struct GMT_CTRL *GMT, char *pSrcSRS, char *pDstSRS, int n_pts,
75 double *xi, double *yi, double *zi, bool insitu, double *xo, double *yo, double *zo) {
76 /* pSrcSRS and pDstSRS are pointers to strings defining the Source and Destination
77 Referencing System. The SRS can be a +proj Proj.4 string, a WKT, a EPSG:n code or a filename with a WKT (?).
78 n_pts is the number of points to be transformed.
79 xi,yi,zi are pointers to arrays of n_pts points. If only 2D transform is wanted, pass zi = NULL.
80 insitu, is a boolean stating if the transformed points will overwrite the input data in xi,yi,zi (true)
81 or, when false, the output is stored in xo,yo,zo. In this later case, it's user responsibility
82 to allocate the xo,yo[,zo] arrays with same size as xi,yi[,zi].
83 */
84 OGRCoordinateTransformationH hCT = gmt_OGRCoordinateTransformation(GMT, pSrcSRS, pDstSRS);
85
86 if (insitu)
87 OCTTransform(hCT, n_pts, xi, yi, zi);
88 else {
89 int n;
90 for (n = 0; n < n_pts; n++) {
91 xo[n] = xi[n];
92 yo[n] = yi[n];
93 }
94 if (zi) { /* If z component too */
95 for (n = 0; n < n_pts; n++) zo[n] = zi[n];
96 }
97 OCTTransform(hCT, n_pts, xo, yo, zo);
98 }
99
100 OCTDestroyCoordinateTransformation(hCT);
101 return (GMT_NOERROR);
102 }
103
gmt_ogrproj_one_pt(OGRCoordinateTransformationH hCT,double * xi,double * yi,double * zi)104 void gmt_ogrproj_one_pt(OGRCoordinateTransformationH hCT, double *xi, double *yi, double *zi) {
105 /* Suitable to call on a rec-by-rec basis because *hCT must have been initiated outside.
106 Again, *zi may be NULL for the 2D case and after last point the hCT GDAL object must be cleaned by a call to
107 OCTDestroyCoordinateTransformation(hCT);
108 */
109 OCTTransform(hCT, 1, xi, yi, zi);
110 }
111
gmt_proj4_fwd(struct GMT_CTRL * GMT,double xi,double yi,double * xo,double * yo)112 void gmt_proj4_fwd(struct GMT_CTRL *GMT, double xi, double yi, double *xo, double *yo) {
113 /* Function that have the same signature as the GMT coordinate transforms */
114 *xo = xi; *yo = yi;
115 gmt_ogrproj_one_pt(GMT->current.gdal_read_in.hCT_fwd, xo, yo, NULL);
116 }
117
gmt_proj4_inv(struct GMT_CTRL * GMT,double * xi,double * yi,double xo,double yo)118 void gmt_proj4_inv(struct GMT_CTRL *GMT, double *xi, double *yi, double xo, double yo) {
119 /* The inverse transform */
120 *xi = xo; *yi = yo;
121 gmt_ogrproj_one_pt(GMT->current.gdal_read_in.hCT_inv, xi, yi, NULL);
122 }
123