Using Pyproj to Transform Geographical Projection Coordinate System

The reference ellipsoids of the two coordinate systems are different, and the values of different coordinate systems of a point on the ground are different. The coordinate systems used by different departments are often inconsistent, so they can only be used after conversion. For example, the location of the Beijing observation site currently in use is based on GPS. The geographic coordinate system used by GPS is GCS_WGS_1984, so the geographic coordinate system of its coordinates is also GCS_WGS_1984. If you need to display these points on the map on the Web, Web The projected coordinate system WGS_1984_Web_Mercator_Auxiliary_Sphere at the end requires the operation of converting the geographic coordinate system to the projected coordinate system.

Regarding the relationship between the geographic coordinate system and the projected coordinate system, I won’t go into details here, but I will introduce WKID. The full English name of WKID is Well Known ID, which is a well-known number. This number is something that everyone sits down to discuss, agree and agree with, and it is unique in detail. Many coordinate systems have their own WKID, just like everyone has their own ID number, which is determined from birth. Even if the name is changed, it can still be determined by the ID number. This is the use of spatial data, Conversion, sharing, etc. play a key role.

For example, the following table is a specific example of the WKID format of GCS_WGS_1984 and the projection file of a certain point:

The projected coordinates that need to be converted also have its WKID. For example, the following is a specific example of WGS_1984_Web_Mercator_Auxiliary_Sphere:

If you don’t know the WKID of the coordinate system you need to convert, you can query from the following two websites:

Geographical coordinate system WKID: https://developers.arcgis.com/javascript/3/jshelp/gcs.htm

Projection coordinate system WKID: https://developers.arcgis.com/javascript/3/jshelp/pcs.htm

Here are some of the WKIDs of geographic coordinate systems commonly used in my country:

There are many tools for coordinate system conversion, among which the more commonly used are Proj.4 related libraries. If we use Python for coordinate conversion, there are advanced Pyproj third-party libraries that can be used. The document address is as follows:
http://jswhit.github.io/pyproj/

This library is very simple, we only need to master one of the main functions:

transform(p1, p2, x, y, z=None, radians=False)

Example: x2, y2, z2 =transform(p1, p2, x1, y1, z1, radians=False)

This function represents the coordinate conversion between the p1 coordinate system and the p2 coordinate system. x1, y1, z1 are the coordinates defined by the p1 coordinate system, and z is the height in meters. X2, y2, z3 are the coordinates defined by the p2 coordinate system, which are returned after conversion, and z1=none by default. The Radians parameter indicates whether to return the value in radians.

Let's perform the coordinate conversion of the observation station in Beijing.

# Projection transformation

def proj_trans():
     # Read latitude and longitude
     data = pd.read_excel(u"D:/Visualization/python/file/location.xlsx")

     lon = data.lon.values
     lat = data.lat.values
     print lon, lat

     p1 = pyproj.Proj(init="epsg:4326") # Define the data geographic coordinate system
     p2 = pyproj.Proj(init="epsg:3857") # Define conversion projection coordinate system
     x1, y1 = p1(lon, lat)
     x2, y2 = pyproj.transform(p1, p2, x1, y1, radians=True)
         print x2, y2

In the above code, you need to pay attention to the need to define the geographic coordinate system of the data and the projected coordinate system after the conversion before the conversion, so that the conversion can be carried out with a purpose.

The results before and after conversion are as follows:

Let's use the Project tool in ArcMap to experiment and compare the results of the experiment. The following is the result of the projection conversion in ArcMap, which is basically the same as the above result.