Go to a post about gdal Coordinate Transformation

Source post: http://www.gissky.net/Article/647.htm

1. Coordinate Conversion
Coordinate conversion is always a headache. I don't know if you have written the program code for coordinate conversion. If you have written it, you will understand it. That long string of formulas and the braces that appear in the group are terrible! Only when you have written such a thing will you respect the home we live in. Only then will we find that the Earth is so complex and so great, and the mathematics and geographer of mankind are so great and so great. Fear will make you afraid to litter and spam in the future ~~ Oh, that's big. I am afraid to touch those symbols and formulas. However, since we step into the GIS threshold, we have to face the abyss of coordinates. Forget it, find a wheel that has already been made by a family to fight for your own car. As for how many steel wires and materials the wheel has, how can we make it relatively round, I don't have to worry about it anymore.
There is a kind of brand wheel called proj that I have never heard. Anyway, it is almost the best thing about Coordinate Transformation in the Open Source world. For C/C ++ people, we are lucky because this library can be called directly. If you are a Java person on the Linux platform, you are lucky to use JNI for indirect calls. If you are a Java person on Windows, it may be unfortunate. After two days of hard work, my boss and I did not compile successfully. I checked countless posts and found no solutions. (If you have succeeded, please let me know ~~~). For amateurs, I thought we were unfortunate the last two days, but today I think God prefers python, so I found this and this. Oh ~~~
For pyproj, I am also "unfamiliar" :-). Then let's see how the Coordinate Conversion in OSR works. Translate the official tutorial first. If the translation is incorrect, please look back at the original English version.

1.1. Official tutorials

1.1.1. Overview
The ogrspatialreference and ogrcoordinatetransformation classes provide services for describing Coordinate Systems (projections and baselines) and the conversion between coordinate systems. These services have models and WKT descriptions in the OpenGIS coordinate conversion instructions. Some background information about the OpenGIS coordinate system and service can be found in the simple element (Simple Features for com) of COM and the spatial reference system abstract model document (which can be obtained from opengis.org. Geotiff projection transformation list can also be used to understand the projection formula in WKT. The epsg test and learning web page is also a useful resource.

1.1.2. Define a geographic coordinate system
The geographic coordinate system is encapsulated into the ogrspatialreference class. There are several ways to initialize the ogrspatialreference object to form a legal coordinate system. There are two major coordinate systems. One is geographical coordinates (expressed by longitude and latitude ). The other is projection coordinates (such as UTM, which are measured in meters and other units ).
A geographic coordinate contains a datum (it contains the elliptical body and the anti-flat rate described by the long half axis), a primary Meridian (generally Greenwich Mean), and an angle unit (generally degree ). The following code initializes a geographic coordinate system. It provides all information about the geographical coordinate system of user-defined names.

Toggle line numbers toggle line numbers
1 ogrspatialreference osrs;
2 osrs. setgeogcs ("My geographic coordinate system ",
3 "wgs_1984 ",
4 "My WGS84 spheroid ",
5 srs_wgs84_semimajor, srs_wgs84_invflattening,
6 "Greenwich", 0.0,
7 "degree", srs_ua_degree_conv );
Among these values, "My geographic coordinate system", "My WGS84 spheroid", "Greenwich", and "degree" are not keywords. But it is used for display to users. However, the ing name "wgs_1984" is often used to define the ing. What values can be used and what is not used are well-regulated. It cannot be written in disorder.
Ogrspatialreference already supports some standard coordinate systems. For example, "nad27", "nad83", "wgs72" and "WGS84 ". To build them, you only need to use a function setwellknowngeogcs ().

--------------------------------------------------------------------------------

Osrs. setwellknowngeogcs ("WGS84 ");

--------------------------------------------------------------------------------

If an epsg database exists, the geographic coordinate system of all epsg databases can be represented by GCS.

--------------------------------------------------------------------------------

Osrs. setwellknowngeogcs ("epsg: 4326 ");

--------------------------------------------------------------------------------

In the expressions of Coordinate Systems, WKT is a bond, through which various expressions can be exchanged. An ogrspatialreference can be initialized with a WKT or converted to a WKT expression.

Toggle line numbers toggle line numbers
1 char * pszwkt = NULL;
2 osrs. setwellknowngeogcs ("WGS84 ");
3 osrs. exporttowkt (& pszwkt );
4 printf ("% s/n", pszwkt );
For example, a WKT:

Toggle line numbers toggle line numbers
1 geogcs ["WGS 84", datum ["wgs_1984", spheroid ["WGS 84", 6378137,298.257223563,
2 Authority ["epsg", 7030], towgs84 [6326,], authority ["epsg",],
3 primem ["Greenwich", 0, authority ["epsg", 8901], Unit ["dmsh", 0.0174532925199433,
4 authority ["epsg", 9108], axis ["Lat", North], axis ["long", East], authority ["epsg ",
5 4326]
Or let it look better:

Toggle line numbers toggle line numbers
1 geogcs ["WGS 84 ",
2 datum ["wgs_1984 ",
3 spheroid ["WGS 84", 6378137,298.257223563,
4 authority ["epsg", 7030],
5 towgs84 [0, 0, 0, 0, 0],
6 authority ["epsg", 6326],
7 primem ["Greenwich", 0, authority ["epsg", 8901],
8 Unit ["dmsh", 0.0174532925199433, authority ["epsg", 9108],
9 axis ["Lat", North],
10 axis ["long", East],
11 authority ["epsg", 4326]
The ogrspatialreference: importfromwkt () method can be used to set the WKT coordinate system to ogrspatialreference.

 

1.1.3. Define a projection system
A projection coordinate system (such as UTM) is defined based on a geographic coordinate system to facilitate conversion between linear units and angle positions. The following code defines a utm17 projection coordinate system based on the geographic coordinate system under the WGS84.

Toggle line numbers toggle line numbers
1 ogrspatialreference osrs;
2 osrs. setprojcs ("UTM 17 (WGS84) in northern hemisphere .");
3 osrs. setwellknowngeogcs ("WGS84 ");
4 osrs. setutm (17, true );
Call setprojcs () to set up a projection coordinate system with a user-defined name and determine that the system has been projected. Setwellknowngeogcs () allocates a geographic coordinate system, and setutm () sets the details of the projection conversion parameters. A valid definition cannot be created until the preceding steps are completed. However, this object will automatically orchestrate the internal expression when needed to maintain legitimacy.
Pay special attention to the step sequence when defining objects !!
The above definition defines a WKT as below. Note that utm117 is extended to the horizontal axis mocato to define the UTM parameter.

Toggle line numbers toggle line numbers
1 projcs ["UTM 17 (WGS84) in Northern hemisloud .",
2 geogcs ["WGS 84 ",
3 datum ["wgs_1984 ",
4 spheroid ["WGS 84", 6378137,298.257223563,
5 Authority ["epsg", 7030],
6 towgs84 [0, 0, 0, 0, 0],
7 Authority ["epsg", 6326],
8 primem ["Greenwich", 0, authority ["epsg", 8901],
9 Unit ["dmsh", 0.0174532925199433, authority ["epsg", 9108],
10 axis ["Lat", North],
11 axis ["long", East],
12 authority ["epsg", 4326],
13 projection ["transverse_mercator"],
14 parameter ["latitude_of_origin", 0],
15 parameter ["central_meridian",-81],
16 parameter ["scale_factor", 0.9996],
17 parameter ["false_easting", 500000],
18 parameter ["false_northing", 0]
Different projections have different functions to set parameters (settm () (Transverse Mercator), setsvc () (Lambert conformal CONIC), and setmercator ())

1.1.4. query the Coordinate System
Once ogrspatialreference is created, multiple types of information can be queried. Use ogrspatialreference: isprojected () and ogrspatialreference: isgeographic () to check whether it is a projection coordinate or a geographical coordinate. Ogrspatialreference: getsemimajor (), ogrspatialreference: getsemiminor () and ogrspatialreference: getinvflattening () can be used to obtain information about the elliptical body. Ogrspatialreference: getattrvalue () can be used to obtain the name expression strings of projcs, geogcs, datum, spheroid, and projection. Ogrspatialreference: getprojparm () can be used to obtain projection parameters. Ogrspatialreference: getlinearunits () can be used to obtain linear units and convert them to meters.
The following code demonstrates the use of getattrvalue () to obtain projection and the use of getprojparm () getattrvalue. The defined projection parameters (such as srs_pp_central_meridian) should be used when getprojparm () is used to obtain the projection parameters. The code for setting different projection methods in ogrspatialreference. cpp can be used to find which projection is used and which parameter.

Toggle line numbers toggle line numbers
1 const char * pszprojection = posrs-> getattrvalue ("projection ");
2 If (pszprojection = NULL)
3 {
4 If (posrs-> isgeographic ())
5 sprintf (szproj4 + strlen (szproj4), "+ proj = longlat ");
6 else
7 sprintf (szproj4 + strlen (szproj4), "unknown ");
8}
9 else if (equal (pszprojection, srs_pt_cylindrical_assist_area ))
10 {
11 sprintf (szproj4 + strlen (szproj4 ),
12 "+ proj = CEA + lon_0 = %. 9f + lat_ts = %. 9f + x_0 = %. 3f + y_0 = %. 3f ",
13 posrs-> getprojparm (srs_pp_central_meridian, 0.0 ),
14 posrs-> getprojparm (srs_pp_standard_parallel_1, 0.0 ),
15 posrs-> getprojparm (srs_pp_false_easting, 0.0 ),
16 posrs-> getprojparm (srs_pp_false_northing, 0.0 ));
17}
18...

1.1.5. Coordinate Transformation
The ogrcoordinatetransformation class is used to convert coordinates between different coordinate systems. The new conversion object is created by ogrcoordinatetransformation. After creation, the ogrcoordinatetransformation: Transform () method can be used to convert points in different coordinate systems.

Toggle line numbers toggle line numbers
1 ogrspatialreference osourcesrs, otargetsrs;
2 ogrcoordinatetransformation * poct;
3 double X, Y;
4 osourcesrs. importfromepsg (atoi (papszargv [I + 1]);
5 otargetsrs. importfromepsg (atoi (papszargv [I + 2]);
6 poct = ogrcreatecoordinatetransformation (& osourcesrs,
7 & otargetsrs );
8 x = atof (papszargv [I + 3]);
9 y = atof (papszargv [I + 4]);
10 if (poct = NULL |! Poct-> transform (1, & X, & Y ))
11 printf ("transformation failed./N ");
12 else
13 printf ("(% F, % F)-> (% F, % F)/n ",
14 atof (papszargv [I + 3]),
15 atof (papszargv [I + 4]),
16 x, y );
There are two causes for an error during conversion. First, ogrcreatecoordinatetransformation () may fail. Generally, the program considers that the two specified coordinate systems cannot establish a conversion relationship. This may be because the projection proj is temporarily not supported, so that there is no known relationship between the two different baselines. Or a projection is not properly defined. If ogrcreatecoordinatetransformation () fails, null is returned.
The ogrcoordinatetransformation: Transform () method may also fail. This may be due to a problem formed after the above problems are accumulated. It is also the consequence of "omitting an operation because one or more vertices are not defined in the value. If the transform () function succeeds, true is returned. If a point fails to be converted, the remaining points are maintained in an incorrect state.
The coordinate conversion service can actually process 3D points. In addition, the Service adjusts the Elevation Based on the elevation difference between the elliptical body and the datum plane. In the future, there may also be applications for point-to-point conversion on different vector baselines. If there is no Z value, it is assumed that all vertices are on the geolevel.
The following example shows how to easily create a geographic coordinate system to a projection Coordinate System Based on the same geographic coordinate system. And convert between two coordinate systems.

Toggle line numbers toggle line numbers
1 ogrspatialreference outm, * polatlong;
2 ogrcoordinatetransformation * potransform;
3 outm. setprojcs ("UTM 17/WGS84 ");
4 outm. setwellknowngeogcs ("WGS84 ");
5 outm. setutm (17 );
6 polatlong = outm. clonegeogcs ();
7 potransform = ogrcreatecoordinatetransformation (& outm, polatlong );
8 If (potransform = NULL)
9 {
10...
11}
12...
13 if (! Potransform-> transform (npoints, x, y, z ))
14...

1.1.6. API binding
C interfaces are defined in ogr_srs_api.h. python is bound to the OSR. py module. Some C ++ methods are not defined in C and Python bindings.
Python bindings (C binding skipped)

Toggle line numbers toggle line numbers
1 class OSR. spatialreference
2 def _ init _ (self, OBJ = none ):
3 def importfromwkt (self, WKT ):
4 def exporttowkt (Self ):
5 def importfromepsg (self, code ):
6 def isgeographic (Self ):
7 def isprojected (Self ):
8 def getattrvalue (self, name, child = 0 ):
9 def setattrvalue (self, name, value ):
10 def setwellknowngeogcs (self, name ):
11 def setprojcs (self, name = "unnamed "):
12 def issamegeogcs (self, other ):
13 def issame (self, other ):
14 def setlinearunits (self, units_name, to_meters ):
15 def setutm (self, zone, is_north = 1 ):
16 class coordinatetransformation:
17 def _ init _ (self, source, target ):
18 def transformpoint (self, x, y, z = 0 ):
19 def transformpoints (self, points ):

1.1.7. Interface implementation
The implementation of the ogrcoordinatetransformation service is built on the proj4 library, which was first written by Gerald evenden of USGS.

 

1.2. My Experiments

1.2.1. Create and compare

Toggle line numbers toggle line numbers
1 >>> import gdal
2 >>> import OSR
3 >>> dataset = gdal. Open ("E:/gisdata/gtif/spot. tif ")
Obtain the spatial reference from the dataset and create a spatialreference object.

Toggle line numbers toggle line numbers
1 >>> sr = dataset. getprojectionref ()
2 >>> SR
3 'projcs ["unnamed", geogcs ["nad27", datum ["north_american_datum_1927", spheroid ["linoleic"
4 rke 1866 ", 6378206.4, 294.9786982139006, authority [" epsg "," 7008 "], authority [" epsg"
5, "6267"], primem ["Greenwich", 0], Unit ["degree", 0.0174532925199433], authority ["EPS
6g "," 4267 "], projection [" transverse_mercator "], parameter [" latitude_of_origin ", 0],
7 parameter ["central_meridian",-105], parameter ["scale_factor", 0.9996], parameter ["F
8 alse_easting ", 500000], parameter [" false_northing ", 0], Unit [" metre ", 1, authority [" EP
9 SG "," 9001 "], authority [" epsg "," 26713 "]'
10 >>> osrobj = OSR. spatialreference ()
11 >>> osrobj. importfromwkt (SR)
12 0
13 >>> osrobj. exporttowkt ()
14 'projcs ["unnamed", geogcs ["nad27", datum ["north_american_datum_1927", spheroid ["linoleic"
15 rke 1866 ", 6378206.4, 294.9786982139006, authority [" epsg "," 7008 "], authority [" epsg"
16, "6267"], primem ["Greenwich", 0], Unit ["degree", 0.0174532925199433], authority ["EPS
17G "," 4267 "], projection [" transverse_mercator "], parameter [" latitude_of_origin ", 0],
18 parameter ["central_meridian",-105], parameter ["scale_factor", 0.9996], parameter ["F
19 alse_easting ", 500000], parameter [" false_northing ", 0], Unit [" metre ", 1, authority [" EP
20 SG "," 9001 "], authority [" epsg "," 26713 "]'
21 >>> osrobj. isgeographic ()
22 0
23 >>> osrobj. isprojected ()
24 1
Error: EOF in multi-line statement

Obtain another one for comparison.

Toggle line numbers toggle line numbers
1 >>> dataset2 = gdal. Open ("E:/gisdata/gtif/uparea. tif ")
2 >>> Sr2 = dataset2.getprojectionref ()
3 >>> Sr2
4 'projcs ["unnamed", geogcs ["nad27", datum ["north_american_datum_1927", spheroid ["linoleic"
5 rke 1866 ", 6378206.4, 294.9786982139006, authority [" epsg "," 7008 "], authority [" epsg"
6, "6267"], primem ["Greenwich", 0], Unit ["degree", 0.0174532925199433], authority ["EPS
7G "," 4267 "], projection [" transverse_mercator "], parameter [" latitude_of_origin ", 0],
8 parameter ["central_meridian",-105], parameter ["scale_factor", 0.9996], parameter ["F
9 alse_easting ", 500000], parameter [" false_northing ", 0], Unit [" metre ", 1, authority [" EP
10 SG "," 9001 "], authority [" epsg "," 26713 "]'
11 >>> osrobj2 = OSR. spatialreference ()
12 >>> osrobj2.importfromwkt (Sr2)
13 0
14 >>> osrobj2.issame (osrobj)
15 1
Error: EOF in multi-line statement

Create a geographic coordinate system and compare it with an existing Coordinate System

Toggle line numbers toggle line numbers
1 >>> osrobj3 = OSR. spatialreference ()
2 >>> osrobj3.setwellknowngeogcs ("WGS84 ")
3 0
4 >>> osrobj3.issame (osrobj2)
5 0
6 >>> osrobj3.issame (osrobj)
7 0
8 >>> osrobj3.exporttowkt ()
9 'geogcs ["WGS 84", datum ["wgs_1984", spheroid ["WGS 84", 6378137,298.257223563, author
10 ity ["epsg", "7030"], towgs84 [6326, 0], authority ["epsg", ""], primem ["GRE
11 enwich ", 0, authority [" epsg "," 8901 "], Unit [" degree ", 0.0174532925199433, authority ["
12 epsg "," 9108 "], axis [" Lat ", North], axis [" long ", East], authority [" epsg "," 4326 "]'
13 >>> osrobj3.isgeographic ()
14 1

1.2.2. Coordinate Transformation between geographic coordinate system and projection Coordinate System
The coordinate range of the SPOT image projection is as follows:

Toggle line numbers toggle line numbers
1 In projected or local coordinates
2 left: 590000.000000
3 right: 609000.000000
4 top: 4928000.000000
5 bottom: 4914000.000000
Result of coordinatetransformation conversion:

Toggle line numbers toggle line numbers
1 >>> Ct = OSR. coordinatetransformation (osrobj, osrobj3)
2 >>> CT. transformpoint (5920.,00008000)
3 (-103.86789483706976, 44.501658895816213, 2.1979236960925293e-007)
4 >>> CT. transformpoint (609000,00008000)
5 (-103.6289570318164, 44.499040134967487, 2.2072345018386841e-007)
6 >>> CT. transformpoint (609000,4914000)
7 (-103.63190060724861, 44.373035580280714, 2.2072345018386841e-007)
8 >>> CT. transformpoint (59213,4914000)
9 (-103.87032571673741, 44.375642923223602, 2.1979236960925293e-007)
10 >>>
The geographic coordinate range in ArcGIS is as follows:

Toggle line numbers toggle line numbers
1 in decimal degrees
2 West:-103.870326
3 East:-103.628957
4 north: 44.501659
5 South: 44.373036
There are some differences between the products calculated using proj and ArcGIS (one is accurate and the other is 0.01 degrees different ). However, I think this error should be allowed.