.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/projections.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_projections.py: Geographic Coordinates ====================== Most gridders and processing methods in Verde operate under the assumption that the data coordinates are Cartesian. To process data in geographic (longitude and latitude) coordinates, we must first project them. There are different ways of doing this in Python but most of them rely on the `PROJ library `__. We'll use `pyproj `__ to access PROJ directly and handle the projection operations. .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: default import cartopy.crs as ccrs import matplotlib.pyplot as plt import numpy as np import pyproj import verde as vd .. GENERATED FROM PYTHON SOURCE LINES 27-30 With pyproj, we can create functions that will project our coordinates to and from different coordinate systems. For our Baja California bathymetry data, we'll use a Mercator projection. .. GENERATED FROM PYTHON SOURCE LINES 30-35 .. code-block:: default data = vd.datasets.fetch_baja_bathymetry() # We're choosing the latitude of true scale as the mean latitude of our dataset projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) .. GENERATED FROM PYTHON SOURCE LINES 36-38 The Proj object is a callable (meaning that it behaves like a function) that will take longitude and latitude and return easting and northing coordinates. .. GENERATED FROM PYTHON SOURCE LINES 38-43 .. code-block:: default # pyproj doesn't play well with Pandas so we need to convert to numpy arrays proj_coords = projection(data.longitude.values, data.latitude.values) print(proj_coords) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none (array([-11749316.65303046, -11748999.90777981, -11748682.14077029, ..., -11749178.7155826 , -11749531.22239381, -11749882.70744613]), array([2905766.05183735, 2905457.83817751, 2905148.48639437, ..., 2579612.89349502, 2580017.48641078, 2580422.09116873])) .. GENERATED FROM PYTHON SOURCE LINES 44-45 We can plot our projected coordinates using matplotlib. .. GENERATED FROM PYTHON SOURCE LINES 45-56 .. code-block:: default plt.figure(figsize=(7, 6)) plt.title("Projected coordinates of bathymetry measurements") # Plot the bathymetry data locations as black dots plt.plot(proj_coords[0], proj_coords[1], ".k", markersize=0.5) plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal") plt.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_projections_001.png :alt: Projected coordinates of bathymetry measurements :srcset: /tutorials/images/sphx_glr_projections_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 57-63 Cartesian grids --------------- Now we can use :class:`verde.BlockReduce` and :class:`verde.Spline` on our projected coordinates. We'll specify the desired grid spacing as degrees and convert it to Cartesian using the 1 degree approx. 111 km rule-of-thumb. .. GENERATED FROM PYTHON SOURCE LINES 63-68 .. code-block:: default spacing = 10 / 60 reducer = vd.BlockReduce(np.median, spacing=spacing * 111e3) filter_coords, filter_bathy = reducer.filter(proj_coords, data.bathymetry_m) spline = vd.Spline().fit(filter_coords, filter_bathy) .. GENERATED FROM PYTHON SOURCE LINES 69-71 If we now call :meth:`verde.Spline.grid` we'll get back a grid evenly spaced in projected Cartesian coordinates. .. GENERATED FROM PYTHON SOURCE LINES 71-75 .. code-block:: default grid = spline.grid(spacing=spacing * 111e3, data_names="bathymetry") print("Cartesian grid:") print(grid) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Cartesian grid: Dimensions: (northing: 61, easting: 54) Coordinates: * easting (easting) float64 -1.175e+07 -1.173e+07 ... -1.077e+07 * northing (northing) float64 2.074e+06 2.093e+06 ... 3.171e+06 3.19e+06 Data variables: bathymetry (northing, easting) float64 -3.635e+03 -3.727e+03 ... 8.87e+03 Attributes: metadata: Generated by Spline(mindist=0) .. GENERATED FROM PYTHON SOURCE LINES 76-78 We'll mask our grid using :func:`verde.distance_mask` to get rid of all the spurious solutions far away from the data points. .. GENERATED FROM PYTHON SOURCE LINES 78-92 .. code-block:: default grid = vd.distance_mask(proj_coords, maxdist=30e3, grid=grid) plt.figure(figsize=(7, 6)) plt.title("Gridded bathymetry in Cartesian coordinates") pc = grid.bathymetry.plot.pcolormesh(cmap="viridis", vmax=0, add_colorbar=False) plt.colorbar(pc).set_label("bathymetry (m)") plt.plot(filter_coords[0], filter_coords[1], ".k", markersize=0.5) plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal") plt.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_projections_002.png :alt: projections :srcset: /tutorials/images/sphx_glr_projections_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 93-106 Geographic grids ---------------- The Cartesian grid that we generated won't be evenly spaced if we convert the coordinates back to geographic latitude and longitude. Verde gridders allow you to generate an evenly spaced grid in geographic coordinates through the ``projection`` argument of the :meth:`~verde.base.BaseGridder.grid` method. By providing a projection function (like our pyproj ``projection`` object), Verde will generate coordinates for a regular grid and then pass them through the projection function before predicting data values. This way, you can generate a grid in a coordinate system other than the one you used to fit the spline. .. GENERATED FROM PYTHON SOURCE LINES 106-122 .. code-block:: default # Get the geographic bounding region of the data region = vd.get_region((data.longitude, data.latitude)) print("Data region in degrees:", region) # Specify the region and spacing in degrees and a projection function grid_geo = spline.grid( region=region, spacing=spacing, projection=projection, dims=["latitude", "longitude"], data_names="bathymetry", ) print("Geographic grid:") print(grid_geo) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Data region in degrees: (245.0, 254.705, 20.0, 29.99131) Geographic grid: Dimensions: (latitude: 61, longitude: 59) Coordinates: * longitude (longitude) float64 245.0 245.2 245.3 ... 254.4 254.5 254.7 * latitude (latitude) float64 20.0 20.17 20.33 20.5 ... 29.66 29.82 29.99 Data variables: bathymetry (latitude, longitude) float64 -3.621e+03 ... 9.001e+03 Attributes: metadata: Generated by Spline(mindist=0) .. GENERATED FROM PYTHON SOURCE LINES 123-128 Notice that grid has longitude and latitude coordinates and slightly different number of points than the Cartesian grid. The :func:`verde.distance_mask` function also supports the ``projection`` argument and will project the coordinates before calculating distances. .. GENERATED FROM PYTHON SOURCE LINES 128-133 .. code-block:: default grid_geo = vd.distance_mask( (data.longitude, data.latitude), maxdist=30e3, grid=grid_geo, projection=projection ) .. GENERATED FROM PYTHON SOURCE LINES 134-135 Now we can use the Cartopy library to plot our geographic grid. .. GENERATED FROM PYTHON SOURCE LINES 135-146 .. code-block:: default plt.figure(figsize=(7, 6)) ax = plt.axes(projection=ccrs.Mercator()) ax.set_title("Geographic grid of bathymetry") pc = grid_geo.bathymetry.plot.pcolormesh( ax=ax, transform=ccrs.PlateCarree(), vmax=0, zorder=-1, add_colorbar=False ) plt.colorbar(pc).set_label("meters") vd.datasets.setup_baja_bathymetry_map(ax, land=None) plt.show() .. image-sg:: /tutorials/images/sphx_glr_projections_003.png :alt: projections :srcset: /tutorials/images/sphx_glr_projections_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none /home/runner/work/verde/verde/doc/tutorials_src/projections.py:143: UserWarning: All kwargs are being ignored. They are accepted to guarantee backward compatibility. vd.datasets.setup_baja_bathymetry_map(ax, land=None) .. GENERATED FROM PYTHON SOURCE LINES 147-165 Profiles -------- For profiles, things are a bit different. The projection is applied to the input points before coordinates are generated. So the profile will be evenly spaced in *projected coordinates*, not geographic coordinates. This is to avoid issues with calculating distances on a sphere. The coordinates returned by the ``profile`` method will be in geographic coordinates, so projections given to ``profile`` must take an ``inverse`` argument so we can undo the projection. The main utility of using a projection with ``profile`` is being able to pass in points in geographic coordinates and get coordinates back in that same system (making it easier to plot on a map). To generate a profile cutting across our bathymetry data, we can use longitude and latitude points taken from the map above). .. GENERATED FROM PYTHON SOURCE LINES 165-178 .. code-block:: default start = (-114.5, 24.7) end = (-110, 20.5) profile = spline.profile( point1=start, point2=end, size=200, projection=projection, dims=("latitude", "longitude"), data_names=["bathymetry"], ) print(profile) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none latitude longitude distance bathymetry 0 24.700000 -114.500000 0.000000 -4115.540285 1 24.679226 -114.477387 3276.548360 -4397.069051 2 24.658449 -114.454774 6553.096720 -4766.953632 3 24.637668 -114.432161 9829.645080 -5198.947538 4 24.616884 -114.409548 13106.193440 -5639.303390 .. ... ... ... ... 195 20.585669 -110.090452 638926.930221 -3016.467547 196 20.564257 -110.067839 642203.478581 -2969.351779 197 20.542841 -110.045226 645480.026941 -2914.069305 198 20.521422 -110.022613 648756.575301 -2859.557817 199 20.500000 -110.000000 652033.123661 -2818.101948 [200 rows x 4 columns] .. GENERATED FROM PYTHON SOURCE LINES 179-180 Plot the profile location on our geographic grid from above. .. GENERATED FROM PYTHON SOURCE LINES 180-194 .. code-block:: default plt.figure(figsize=(7, 6)) ax = plt.axes(projection=ccrs.Mercator()) ax.set_title("Profile location") pc = grid_geo.bathymetry.plot.pcolormesh( ax=ax, transform=ccrs.PlateCarree(), vmax=0, zorder=-1, add_colorbar=False ) plt.colorbar(pc).set_label("meters") ax.plot(profile.longitude, profile.latitude, "-k", transform=ccrs.PlateCarree()) ax.text(start[0], start[1], "A", transform=ccrs.PlateCarree()) ax.text(end[0], end[1], "B", transform=ccrs.PlateCarree()) vd.datasets.setup_baja_bathymetry_map(ax, land=None) plt.show() .. image-sg:: /tutorials/images/sphx_glr_projections_004.png :alt: projections :srcset: /tutorials/images/sphx_glr_projections_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none /home/runner/work/verde/verde/doc/tutorials_src/projections.py:191: UserWarning: All kwargs are being ignored. They are accepted to guarantee backward compatibility. vd.datasets.setup_baja_bathymetry_map(ax, land=None) .. GENERATED FROM PYTHON SOURCE LINES 195-196 And finally plot the profile. .. GENERATED FROM PYTHON SOURCE LINES 196-207 .. code-block:: default plt.figure(figsize=(8, 3)) ax = plt.axes() ax.set_title("Profile of bathymetry (A-B)") ax.plot(profile.distance, profile.bathymetry, "-k") ax.set_xlabel("Distance (m)") ax.set_ylabel("Bathymetry (m)") ax.set_xlim(profile.distance.min(), profile.distance.max()) ax.grid() plt.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_projections_005.png :alt: Profile of bathymetry (A-B) :srcset: /tutorials/images/sphx_glr_projections_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 2.267 seconds) .. _sphx_glr_download_tutorials_projections.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: projections.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: projections.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_