.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "example_gallery/raster_operations/coordinate_transformations.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_example_gallery_raster_operations_coordinate_transformations.py: Coordinate transformations ========================== .. _example coordinate transformations: Transform a grid from one CRS to another Introduction ------------ A Coordinate Reference System (CRS) is a means of representing a 3D surface (usually of Earth) onto a 2D map. There are many Coordinate Reference Systems, all with their own advantages and discrepencies. There is no one CRS to rule them all, there are only tradeoffs. That said, it is generally recommended to use a locally defined CRS when possible. .. Note :: In gridkit, a CRS is optional. If you represent, say, a generic image as a grid there is no need for a CRS. By default, the CRS property on a grid is None. It can either be supplied when the grid is initiated or by setting the propery on the grid after creation, e.g. ``grid.crs = 4326``. Warped grids ------------ Grids always have straight lines in the CRS they are defined in. When a straigt line is transformed to a different CRS, it often gets warped. To demonstrate this, let's create a grid in WGS84 (epsg code 4326) and transform the coordinates of each cell to CRS UTM N29 (epsg code 32629). .. GENERATED FROM PYTHON SOURCE LINES 23-48 .. code-block:: Python import numpy from gridkit import BoundedRectGrid # Create a new grid grid_wgs84 = BoundedRectGrid( numpy.arange(25 * 25).reshape(25, 25), bounds=(-40, -30, 10, 20), crs=4236, nodata_value=numpy.nan, ).astype( "float32" ) # use float dtype to be able to represent nans as nodata # Plot the grid in it's native CRS import matplotlib.pyplot as plt plt.imshow(grid_wgs84, extent=grid_wgs84.mpl_extent) plt.title("Original grid in WGS84") plt.xlabel("$^\circ$Lon") plt.ylabel("$^\circ$Lat") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_001.png :alt: Original grid in WGS84 :srcset: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 50-53 This is what the grid looks like it it's native CRS. Now let's transform the grid's cell centers .. GENERATED FROM PYTHON SOURCE LINES 54-72 .. code-block:: Python from pyproj.transformer import CRS, Transformer utm_epsg = 32629 transformer = Transformer.from_crs( grid_wgs84.crs, CRS.from_user_input(utm_epsg), always_xy=True ) # UTM zone 28N wgs84_centroids_utm = transformer.transform( *grid_wgs84.centroid().T ) # transform each cell # Plot the cell centers as a scatter plot, color them using their data value plt.scatter(*wgs84_centroids_utm, s=20, c=grid_wgs84.data.ravel(), cmap="viridis") plt.title("WGS84 grid cells transformed to UTM") plt.xlabel("x [metre]") plt.ylabel("y [metre]") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_002.png :alt: WGS84 grid cells transformed to UTM :srcset: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 73-80 Notice how the whole grid is warped. There is nothing we can do to combat the warping. In order to work with a grid in a different CRS it needs to be resampled onto a grid that is straight in the new CRS. That is what happens when we call the ``to_crs()`` method on the grid. A new grid is crated that is straigt in the desired CRS. The values of the new grid cells are interpolated from the transformed values of the source grid. To demonstrate this principle, let's plot the straight grid as red dots on top of the warped one, containing the original data. .. GENERATED FROM PYTHON SOURCE LINES 81-93 .. code-block:: Python grid_utm = grid_wgs84.to_crs(utm_epsg) # resample using to_crs # Plot the cells of the source grid in the new CRS plt.scatter(*wgs84_centroids_utm, s=20, c=grid_wgs84.data.ravel(), cmap="viridis") # Plot the location of the new cells plt.scatter(*grid_utm.centroid().T, s=3, color="red") plt.title("New UTM grid on top of WGS84 grid in UTM") plt.xlabel("x [metre]") plt.ylabel("y [metre]") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_003.png :alt: New UTM grid on top of WGS84 grid in UTM :srcset: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 94-99 Now we showed only the location of the new grid cells, but when ``to_crs()`` was called, the data already got resampled. By default the 'nearest' ``resample_method`` is used. Let's plot the new grid with the transformed WGS84 cell locations on top to see what happened. .. GENERATED FROM PYTHON SOURCE LINES 100-106 .. code-block:: Python plt.imshow(grid_utm, extent=grid_utm.mpl_extent) plt.scatter(*wgs84_centroids_utm, s=3, color="orange") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_004.png :alt: coordinate transformations :srcset: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 107-110 For a smoother result, the 'bilinear' ``resample_method`` can be used. This does lead to artifacts at the boundaries. .. GENERATED FROM PYTHON SOURCE LINES 111-118 .. code-block:: Python grid_utm = grid_wgs84.to_crs(utm_epsg, resample_method="bilinear") plt.imshow(grid_utm, extent=grid_utm.mpl_extent) plt.scatter(*wgs84_centroids_utm, s=3, color="orange") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_005.png :alt: coordinate transformations :srcset: /example_gallery/raster_operations/images/sphx_glr_coordinate_transformations_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-121 ``to_crs()`` is just a shorthand for creating a grid with the same properties as the source grid, but with a different CRS, and calling resample. Any time ``source_grid.resample(target_grid)`` is called, the CRS of the grids are compared. If they are different, the logic described in this example is applied. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.541 seconds) .. _sphx_glr_download_example_gallery_raster_operations_coordinate_transformations.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: coordinate_transformations.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: coordinate_transformations.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: coordinate_transformations.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_