File: build_crs_cf.rst

package info (click to toggle)
python-pyproj 3.7.1-2
  • links: PTS, VCS
  • area: main
  • in suites: sid
  • size: 1,720 kB
  • sloc: python: 13,468; sh: 273; makefile: 90
file content (195 lines) | stat: -rw-r--r-- 5,741 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
.. _build_crs_cf:

Managing CRS to and from CF
============================

http://cfconventions.org/cf-conventions/cf-conventions.html

Exporting CRS to CF
--------------------
When exporting a CRS to the Climate and Forecast (CF) conventions,
you need both the grid mapping as well as the coordinate system.
If you don't use the coordinate system, then you will lose the units
of your projection.


In this example, this is the CRS we will use:

.. code-block:: python

    from pyproj import CRS

    crs = CRS("EPSG:4326")


To get the grid mapping you use :meth:`pyproj.crs.CRS.to_cf`:

.. versionadded:: 2.2.0

.. code-block:: python

    cf_grid_mapping = crs.to_cf()


Contents of `cf_grid_mapping`::


    {'crs_wkt': 'GEOGCRS["WGS 84",DATUM["World Geodetic System '
                ....,ID["EPSG",4326]]',
    'geographic_crs_name': 'WGS 84',
    'grid_mapping_name': 'latitude_longitude',
    'inverse_flattening': 298.257223563,
    'longitude_of_prime_meridian': 0.0,
    'prime_meridian_name': 'Greenwich',
    'reference_ellipsoid_name': 'WGS 84',
    'semi_major_axis': 6378137.0,
    'semi_minor_axis': 6356752.314245179}


To get the coordinate system, you use :meth:`pyproj.crs.CRS.cs_to_cf`:

.. versionadded:: 3.0.0

.. code-block:: python

    cf_coordinate_system = crs.cs_to_cf()


Contents of `cf_coordinate_system`::

    [{'long_name': 'geodetic latitude coordinate',
    'standard_name': 'latitude',
    'units': 'degrees_north',
    'axis': 'Y'},
    {'long_name': 'geodetic longitude coordinate',
    'standard_name': 'longitude',
    'units': 'degrees_east',
    'axis': 'X'}]


Importing CRS from CF
----------------------

When importing a CRS from the Climate and Forecast (CF) conventions,
you need both the grid mapping as well as the coordinate system.
If you don't use the coordinate system, then you will lose the units
of your projection.

.. note:: If the CF `crs_wkt` attribute is available, the coordinate system is
          inside of the WKT and can be used to create the CRS in a single step.

.. warning:: If building from grid mapping, be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224


Build the CRS from CF grid mapping:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In this example, this is the grid mapping and coordinate system we will use::

  variables:
    double x(x) ;
      x:standard_name = "projection_x_coordinate" ;
      x:long_name = "Easting" ;
      x:units = "m" ;
    double y(y) ;
      y:standard_name = "projection_y_coordinate" ;
      y:long_name = "Northing" ;
      y:units = "m" ;
    int crsOSGB ;
      crsOSGB:grid_mapping_name = "transverse_mercator";
      crsOSGB:semi_major_axis = 6377563.396 ;
      crsOSGB:inverse_flattening = 299.3249646 ;
      crsOSGB:longitude_of_prime_meridian = 0.0 ;
      crsOSGB:latitude_of_projection_origin = 49.0 ;
      crsOSGB:longitude_of_central_meridian = -2.0 ;
      crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ;
      crsOSGB:false_easting = 400000.0 ;
      crsOSGB:false_northing = -100000.0 ;

.. note:: If the units are meters as in this example,
          then no further changes are necessary.

.. code-block:: python

    from pyproj import CRS

    crs = CRS.from_cf(
        {
            "grid_mapping_name": "transverse_mercator",
            "semi_major_axis": 6377563.396,
            "inverse_flattening": 299.3249646,
            "longitude_of_prime_meridian": 0.0,
            "latitude_of_projection_origin": 49.0,
            "longitude_of_central_meridian": -2.0,
            "scale_factor_at_central_meridian": 0.9996012717,
            "false_easting": 400000.0,
            "false_northing": -100000.0,
        }
    )


Modify the CRS with coordinate system:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. versionadded:: 3.0.0

.. note:: If the CF `crs_wkt` attribute is available, the coordinate system is
          inside of the WKT and can be used to create the CRS in a single step.

.. warning:: Be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224


In this example, assume everything is the same as above.
However, the units are instead `US_Survey_Foot`::

  variables:
    double x(x) ;
      x:standard_name = "projection_x_coordinate" ;
      x:long_name = "Easting" ;
      x:units = "US_Survey_Foot" ;
    double y(y) ;
      y:standard_name = "projection_y_coordinate" ;
      y:long_name = "Northing" ;
      y:units = "US_Survey_Foot" ;
    ...


In this case, you will need to get the unit conversion factor:

https://github.com/SciTools/cf-units


.. code-block:: python

    from cf_units import Unit
    from pyproj import CRS

    cf_unit = Unit("US_Survey_Foot")
    unit = {
        "type": "LinearUnit",
        "name": "US Survey Foot",
        "conversion_factor": cf_unit.convert(1, "m"),
    }
    cartesian_cs = {
        "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
        "type": "CoordinateSystem",
        "subtype": "Cartesian",
        "axis": [
            {"name": "Easting", "abbreviation": "E", "direction": "east", "unit": unit},
            {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": unit},
        ],
    }
    crs = CRS.from_cf(
        {
            "grid_mapping_name": "transverse_mercator",
            "semi_major_axis": 6377563.396,
            "inverse_flattening": 299.3249646,
            "longitude_of_prime_meridian": 0.0,
            "latitude_of_projection_origin": 49.0,
            "longitude_of_central_meridian": -2.0,
            "scale_factor_at_central_meridian": 0.9996012717,
            "false_easting": 400000.0,
            "false_northing": -100000.0,
        },
        cartesian_cs=cartesian_cs,
    )