Set gridline (if available) as the default grid registration for remo…

…te datasets (#2266)

Co-authored-by: Dongdong Tian <[email protected]>
Co-authored-by: Wei Ji <[email protected]>

pygmt/datasets/earth_age.py

@@ -44,9 +44,7 @@ def load_earth_age(resolution="01d", region=None, registration=None):
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a pixel-registered grid is returned unless only the
-        gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``.
 
     Returns
     -------

pygmt/datasets/earth_free_air_anomaly.py

@@ -44,9 +44,8 @@ def load_earth_free_air_anomaly(resolution="01d", region=None, registration=None
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a pixel-registered grid is returned unless only the
-        gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``
+        for all resolutions except ``"01m"`` which is ``"pixel"`` only.
 
     Returns
     -------

pygmt/datasets/earth_geoid.py

@@ -44,9 +44,7 @@ def load_earth_geoid(resolution="01d", region=None, registration=None):
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a pixel-registered grid is returned unless only the
-        gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``.
 
     Returns
     -------

pygmt/datasets/earth_magnetic_anomaly.py

@@ -49,9 +49,8 @@ def load_earth_magnetic_anomaly(
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a pixel-registered grid is returned unless only the
-        gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``
+        for all resolutions except ``"02m"`` which is ``"pixel"`` only.
 
     mag4km : bool
         Choose the data version to use. The default is ``False``, which is

pygmt/datasets/earth_relief.py

@@ -57,12 +57,8 @@ def load_earth_relief(
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a gridline-registered grid is returned unless only the pixel-registered
-        grid is available.
-
-        **Note**: For GMT 6.3, ``registration=None`` returns a pixel-registered
-        grid by default unless only the gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``
+        for all resolutions except ``"15s"`` which is ``"pixel"`` only.
 
     data_source : str
         Select the source for the Earth relief data.

pygmt/datasets/earth_vertical_gravity_gradient.py

@@ -45,9 +45,8 @@ def load_earth_vertical_gravity_gradient(
 
     registration : str
         Grid registration type. Either ``"pixel"`` for pixel registration or
-        ``"gridline"`` for gridline registration. Default is ``None``, where
-        a pixel-registered grid is returned unless only the
-        gridline-registered grid is available.
+        ``"gridline"`` for gridline registration. Default is ``"gridline"``
+        for all resolutions except ``"01m"`` which is ``"pixel"`` only.
 
     Returns
     -------

pygmt/datasets/load_remote_dataset.py

@@ -230,20 +230,26 @@ def _load_remote_dataset(
     The returned :class:`xarray.DataArray` doesn't support slice operation for
     tiled grids.
     """
-
-    if registration in ("pixel", "gridline", None):
-        # If None, let GMT decide on Pixel/Gridline type
-        reg = f"_{registration[0]}" if registration else ""
+    dataset = datasets[dataset_name]
+    if resolution not in dataset.resolutions.keys():
+        raise GMTInvalidInput(f"Invalid resolution '{resolution}'.")
+    if registration is None:
+        # Check if "gridline" is an available registration for the resolution
+        if "gridline" in dataset.resolutions[resolution].registrations:
+            # Use default of gridline registration if available
+            registration = "gridline"
+        else:
+            registration = "pixel"
+    if registration in ("pixel", "gridline"):
+        reg = f"_{registration[0]}"
     else:
         raise GMTInvalidInput(
             f"Invalid grid registration: '{registration}', should be either "
             "'pixel', 'gridline' or None. Default is None, where a "
-            "pixel-registered grid is returned unless only the "
-            "gridline-registered grid is available."
+            "gridline-registered grid is returned unless only the "
+            "pixel-registered grid is available."
         )
-    dataset = datasets[dataset_name]
-    if resolution not in dataset.resolutions.keys():
-        raise GMTInvalidInput(f"Invalid resolution '{resolution}'.")
+
     if registration and (
         registration not in dataset.resolutions[resolution].registrations
     ):

pygmt/helpers/testing.py

@@ -175,25 +175,26 @@ def download_test_data():
         "@N37W120.earth_relief_03s_g.nc",
         "@N00W090.earth_relief_03m_p.nc",
         "@N00E135.earth_relief_30s_g.nc",
+        "@N00W010.earth_relief_15s_p.nc",  # Specific grid for 15s test
+        "@N04W010.earth_relief_03s_g.nc",  # Specific grid for 03s test
         # Earth synbath relief grids
         "@S15W105.earth_synbath_30s_p.nc",
         # Earth seafloor age grids
         "@earth_age_01d_g",
-        "@S90W180.earth_age_05m_g.nc",  # Specific grid for 05m test
+        "@N00W030.earth_age_01m_g.nc",  # Specific grid for 01m test
         # Earth geoid grids
         "@earth_geoid_01d_g",
-        "@S90W180.earth_geoid_05m_g.nc",  # Specific grid for 05m test
+        "@N00W030.earth_geoid_01m_g.nc",  # Specific grid for 01m test
         # Earth magnetic anomaly grids
         "@earth_mag_01d_g",
-        "@S90W180.earth_mag_05m_g.nc",  # Specific grid for 05m test
+        "@S30W060.earth_mag_02m_p.nc",  # Specific grid for 02m test
         "@earth_mag4km_01d_g",
-        "@S90W180.earth_mag4km_05m_g.nc",  # Specific grid for 05m test
         # Earth free-air anomaly grids
         "@earth_faa_01d_g",
-        "@S90W180.earth_faa_05m_g.nc",  # Specific grid for 05m test
+        "@N00W030.earth_faa_01m_p.nc",  # Specific grid for 01m test
         # Earth vertical gravity gradient grids
         "@earth_vgg_01d_g",
-        "@S90W180.earth_vgg_05m_g.nc",  # Specific grid for 05m test
+        "@N00W030.earth_vgg_01m_p.nc",  # Specific grid for 01m test
         # Other cache files
         "@capitals.gmt",
         "@earth_relief_20m_holes.grd",

pygmt/tests/test_datasets_earth_age.py

@@ -31,12 +31,13 @@ def test_earth_age_01d():
     """
     Test some properties of the earth age 01d data.
     """
-    data = load_earth_age(resolution="01d", registration="gridline")
+    data = load_earth_age(resolution="01d")
     assert data.name == "seafloor_age"
     assert data.attrs["units"] == "Myr"
     assert data.attrs["long_name"] == "age of seafloor crust"
     assert data.attrs["horizontal_datum"] == "WGS84"
     assert data.shape == (181, 361)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
     npt.assert_allclose(data.lon, np.arange(-180, 181, 1))
     npt.assert_allclose(data.min(), 0.167381, rtol=1e-5)
@@ -47,39 +48,21 @@ def test_earth_age_01d_with_region():
     """
     Test loading low-resolution earth age with 'region'.
     """
-    data = load_earth_age(
-        resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
-    )
+    data = load_earth_age(resolution="01d", region=[-10, 10, -5, 5])
     assert data.shape == (11, 21)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
     npt.assert_allclose(data.lon, np.arange(-10, 11, 1))
     npt.assert_allclose(data.min(), 11.293945)
     npt.assert_allclose(data.max(), 125.1189)
 
 
-def test_earth_age_05m_with_region():
-    """
-    Test loading a subregion of high-resolution earth age.
-    """
-    data = load_earth_age(
-        resolution="05m", region=[-50, -40, 20, 30], registration="gridline"
-    )
-    assert data.coords["lat"].data.min() == 20.0
-    assert data.coords["lat"].data.max() == 30.0
-    assert data.coords["lon"].data.min() == -50.0
-    assert data.coords["lon"].data.max() == -40.0
-    npt.assert_allclose(data.data.min(), 0.040000916)
-    npt.assert_allclose(data.data.max(), 46.530003)
-    assert data.sizes["lat"] == 121
-    assert data.sizes["lon"] == 121
-
-
-def test_earth_age_05m_without_region():
+def test_earth_age_01m_without_region():
     """
     Test loading high-resolution earth age without passing 'region'.
     """
     with pytest.raises(GMTInvalidInput):
-        load_earth_age("05m")
+        load_earth_age("01m")
 
 
 def test_earth_age_incorrect_resolution_registration():
@@ -89,3 +72,19 @@ def test_earth_age_incorrect_resolution_registration():
     """
     with pytest.raises(GMTInvalidInput):
         load_earth_age(resolution="01m", region=[0, 1, 3, 5], registration="pixel")
+
+
+def test_earth_age_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has
+    a "gridline" registration.
+    """
+    data = load_earth_age(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (121, 61)
+    assert data.gmt.registration == 0
+    assert data.coords["lat"].data.min() == 3.0
+    assert data.coords["lat"].data.max() == 5.0
+    assert data.coords["lon"].data.min() == -10.0
+    assert data.coords["lon"].data.max() == -9.0
+    npt.assert_allclose(data.min(), 88.63)
+    npt.assert_allclose(data.max(), 125.25)

pygmt/tests/test_datasets_earth_free_air_anomaly.py

@@ -31,12 +31,13 @@ def test_earth_faa_01d():
     """
     Test some properties of the free air anomaly 01d data.
     """
-    data = load_earth_free_air_anomaly(resolution="01d", registration="gridline")
+    data = load_earth_free_air_anomaly(resolution="01d")
     assert data.name == "free_air_anomaly"
     assert data.attrs["long_name"] == "IGPP Global Earth Free-Air Anomaly"
     assert data.attrs["units"] == "mGal"
     assert data.attrs["horizontal_datum"] == "WGS84"
     assert data.shape == (181, 361)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
     npt.assert_allclose(data.lon, np.arange(-180, 181, 1))
     npt.assert_allclose(data.min(), -275.75)
@@ -47,36 +48,35 @@ def test_earth_faa_01d_with_region():
     """
     Test loading low-resolution earth free air anomaly with 'region'.
     """
-    data = load_earth_free_air_anomaly(
-        resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
-    )
+    data = load_earth_free_air_anomaly(resolution="01d", region=[-10, 10, -5, 5])
     assert data.shape == (11, 21)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
     npt.assert_allclose(data.lon, np.arange(-10, 11, 1))
     npt.assert_allclose(data.min(), -58.75)
     npt.assert_allclose(data.max(), 69.524994)
 
 
-def test_earth_faa_05m_with_region():
+def test_earth_faa_01m_without_region():
     """
-    Test loading a subregion of high-resolution earth free air anomaly data.
+    Test loading high-resolution earth free air anomaly without passing
+    'region'.
     """
-    data = load_earth_free_air_anomaly(
-        resolution="05m", region=[-115, -112, 4, 6], registration="gridline"
-    )
-    assert data.shape == (25, 37)
-    assert data.lat.min() == 4
-    assert data.lat.max() == 6
-    assert data.lon.min() == -115
-    assert data.lon.max() == -112
-    npt.assert_allclose(data.min(), -20.5)
-    npt.assert_allclose(data.max(), -3.9500122)
+    with pytest.raises(GMTInvalidInput):
+        load_earth_free_air_anomaly("01m")
 
 
-def test_earth_faa_05m_without_region():
+def test_earth_faa_01m_default_registration():
     """
-    Test loading high-resolution earth free air anomaly without passing
-    'region'.
+    Test that the grid returned by default for the 1 arc-minute resolution has
+    a "pixel" registration.
     """
-    with pytest.raises(GMTInvalidInput):
-        load_earth_free_air_anomaly("05m")
+    data = load_earth_free_air_anomaly(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (120, 60)
+    assert data.gmt.registration == 1
+    npt.assert_allclose(data.coords["lat"].data.min(), 3.008333333)
+    npt.assert_allclose(data.coords["lat"].data.max(), 4.991666666)
+    npt.assert_allclose(data.coords["lon"].data.min(), -9.99166666)
+    npt.assert_allclose(data.coords["lon"].data.max(), -9.00833333)
+    npt.assert_allclose(data.min(), -51)
+    npt.assert_allclose(data.max(), 113.675)

pygmt/tests/test_datasets_earth_geoid.py

@@ -31,12 +31,13 @@ def test_earth_geoid_01d():
     """
     Test some properties of the earth geoid 01d data.
     """
-    data = load_earth_geoid(resolution="01d", registration="gridline")
+    data = load_earth_geoid(resolution="01d")
     assert data.name == "earth_geoid"
     assert data.attrs["units"] == "m"
     assert data.attrs["long_name"] == "EGM2008 Global Earth Geoid"
     assert data.attrs["horizontal_datum"] == "WGS84"
     assert data.shape == (181, 361)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-90, 91, 1))
     npt.assert_allclose(data.lon, np.arange(-180, 181, 1))
     npt.assert_allclose(data.min(), -106.45)
@@ -47,39 +48,21 @@ def test_earth_geoid_01d_with_region():
     """
     Test loading low-resolution earth geoid with 'region'.
     """
-    data = load_earth_geoid(
-        resolution="01d", region=[-10, 10, -5, 5], registration="gridline"
-    )
+    data = load_earth_geoid(resolution="01d", region=[-10, 10, -5, 5])
     assert data.shape == (11, 21)
+    assert data.gmt.registration == 0
     npt.assert_allclose(data.lat, np.arange(-5, 6, 1))
     npt.assert_allclose(data.lon, np.arange(-10, 11, 1))
     npt.assert_allclose(data.min(), 4.87)
     npt.assert_allclose(data.max(), 29.89)
 
 
-def test_earth_geoid_05m_with_region():
-    """
-    Test loading a subregion of high-resolution earth geoid.
-    """
-    data = load_earth_geoid(
-        resolution="05m", region=[-50, -40, 20, 30], registration="gridline"
-    )
-    assert data.coords["lat"].data.min() == 20.0
-    assert data.coords["lat"].data.max() == 30.0
-    assert data.coords["lon"].data.min() == -50.0
-    assert data.coords["lon"].data.max() == -40.0
-    npt.assert_allclose(data.min(), -32.79)
-    npt.assert_allclose(data.max(), 16.57)
-    assert data.sizes["lat"] == 121
-    assert data.sizes["lon"] == 121
-
-
-def test_earth_geoid_05m_without_region():
+def test_earth_geoid_01m_without_region():
     """
     Test loading high-resolution earth geoid without passing 'region'.
     """
     with pytest.raises(GMTInvalidInput):
-        load_earth_geoid("05m")
+        load_earth_geoid("01m")
 
 
 def test_earth_geoid_incorrect_resolution_registration():
@@ -89,3 +72,19 @@ def test_earth_geoid_incorrect_resolution_registration():
     """
     with pytest.raises(GMTInvalidInput):
         load_earth_geoid(resolution="01m", region=[0, 1, 3, 5], registration="pixel")
+
+
+def test_earth_geoid_01m_default_registration():
+    """
+    Test that the grid returned by default for the 1 arc-minute resolution has
+    a "gridline" registration.
+    """
+    data = load_earth_geoid(resolution="01m", region=[-10, -9, 3, 5])
+    assert data.shape == (121, 61)
+    assert data.gmt.registration == 0
+    assert data.coords["lat"].data.min() == 3.0
+    assert data.coords["lat"].data.max() == 5.0
+    assert data.coords["lon"].data.min() == -10.0
+    assert data.coords["lon"].data.max() == -9.0
+    npt.assert_allclose(data.min(), 20.34)
+    npt.assert_allclose(data.max(), 30.039999)