<p>It looks like you simply messed up the order of your min and max of the latitude. And you were looking at the wrong place (-44° instead of +44°).</p>
<pre><code># bathymetry
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

data = Dataset(&quot;ETOPO1_Bed_g_gmt4.grd&quot;,'r')

lon_range = data.variables['x_range'][:]
lat_range = data.variables['y_range'][:]
topo_range = data.variables['z_range'][:]
spacing = data.variables['spacing'][:]
dimension = data.variables['dimension'][:]
z = data.variables['z'][:]
lon_num = dimension[0]
lat_num = dimension[1]

lon = np.linspace(lon_range[0],lon_range[1],dimension[0])
lat = np.linspace(lat_range[0],lat_range[1],dimension[1])

topo = np.reshape(z, (lat_num, lon_num))

# Select region
lon_idx_min = np.abs(lon-10).argmin() # I increased the margin a little
lon_idx_max = np.abs(lon-20).argmin() # to better see where the region is
lat_idx_min = np.abs(lat+48).argmin()
lat_idx_max = np.abs(lat+40).argmin()

plt.imshow(topo[lat_idx_min:lat_idx_max, lon_idx_min:lon_idx_max], vmax=0)
</code></pre>
<p>By using np.linspace you can remove the loop from your code and speed it up significantly.
There is no need to create a meshgrid from lat and lon.</p>

CC BY-SA 4.0
2022-06-15T10:44:16.997Z
greetings fellow meteorologist
Joscha Fregin
CC BY-SA 4.0
2022-06-15T09:38:59.203Z
2022-06-15T09:38:59.203Z
smichel
<p>I am trying to select a region in the Mediterranean Sea from the ETOPO1 data.</p>
<pre><code># bathymetry
from netCDF4 import Dataset

data = Dataset(&quot;ETOPO1_Bed_g_gdal.grd&quot;,'r')
print(data.variables.keys())

# Get data
lon_range = data.variables['x_range'][:]
lat_range = data.variables['y_range'][:]
topo_range = data.variables['z_range'][:]
spacing = data.variables['spacing'][:]
dimension = data.variables['dimension'][:]
z = data.variables['z'][:]
lon_num = dimension[0]
lat_num = dimension[1]

# Prepare array
lon_input = np.zeros(lon_num); lat_input = np.zeros(lat_num)
for i in range(lon_num):
    lon_input[i] = lon_range[0] + i * spacing[0]
for i in range(lat_num):
    lat_input[i] = lat_range[0] + i * spacing[1]

# Create 2D array
lon, lat = np.meshgrid(lon_input, lat_input)

# Convert 2D array from 1D array for z value
topo = np.reshape(z, (lat_num, lon_num))

# Select region
lon_idx_min = np.abs(lon_input-12.6).argmin()
lon_idx_max = np.abs(lon_input-13.4).argmin()
lat_idx_min = np.abs(lat_input-43.8).argmin()
lat_idx_max = np.abs(lat_input-44.2).argmin()

# get data
topo[lon_idx_min:lon_idx_max, lat_idx_min:lat_idx_max].data
</code></pre>
<p>This returns an empty array:</p>
<pre><code>&gt;&gt;array([], shape=(0, 24), dtype=float64)
</code></pre>
<p>One main questions:</p>
<p><strong>1 - What I am doing wrong?</strong>
When I check the lon lat values I am slicing (lat_input[lat_idx_max] etc) they make sense. However, lon_idx_max and lon_idx_min are above len(topo) (11556 and 11604 vs 10801). I guess this is the problem, but I cannot find the cause of it.</p>
<p>One bonus question:</p>
<p>2 - Why if I flip lon and lat I get some data? Note that the data do not make sense (-4000ish meters in the Adriatic Sea...)</p>

CC BY-SA 4.0
2022-06-22T10:46:25.580Z
hi, where does x_range come from, because when i try to rebuild your example python doesn't know it for me
Weiss
CC BY-SA 4.0
2022-06-22T23:28:22.820Z
when I print the variables of the nc file, it is there along y_range and z_range
shamalaia
<oceanography><python><topography>
ETOPO1 region selection (in python)
shamalaia
# ETOPO1 region selection (in python)

(Asked by: username_0 on 2022-06-12)

I am trying to select a region in the Mediterranean Sea from the ETOPO1 data.

```
# bathymetry
from netCDF4 import Dataset

data = Dataset("ETOPO1_Bed_g_gdal.grd",'r')
print(data.variables.keys())

# Get data
lon_range = data.variables['x_range'][:]
lat_range = data.variables['y_range'][:]
topo_range = data.variables['z_range'][:]
spacing = data.variables['spacing'][:]
dimension = data.variables['dimension'][:]
z = data.variables['z'][:]
lon_num = dimension[0]
lat_num = dimension[1]

# Prepare array
lon_input = np.zeros(lon_num); lat_input = np.zeros(lat_num)
for i in range(lon_num):
    lon_input[i] = lon_range[0] + i * spacing[0]
for i in range(lat_num):
    lat_input[i] = lat_range[0] + i * spacing[1]

# Create 2D array
lon, lat = np.meshgrid(lon_input, lat_input)

# Convert 2D array from 1D array for z value
topo = np.reshape(z, (lat_num, lon_num))

# Select region
lon_idx_min = np.abs(lon_input-12.6).argmin()
lon_idx_max = np.abs(lon_input-13.4).argmin()
lat_idx_min = np.abs(lat_input-43.8).argmin()
lat_idx_max = np.abs(lat_input-44.2).argmin()

# get data
topo[lon_idx_min:lon_idx_max, lat_idx_min:lat_idx_max].data

```

This returns an empty array:

```
>>array([], shape=(0, 24), dtype=float64)

```

One main questions:

**1 - What I am doing wrong?**
When I check the lon lat values I am slicing (lat\_input[lat\_idx\_max] etc) they make sense. However, lon\_idx\_max and lon\_idx\_min are above len(topo) (11556 and 11604 vs 10801). I guess this is the problem, but I cannot find the cause of it.

One bonus question:

2 - Why if I flip lon and lat I get some data? Note that the data do not make sense (-4000ish meters in the Adriatic Sea...)

**Tags:** oceanography, python, topography

> **Comments:**
> **username_1**: hi, where does x_range come from, because when i try to rebuild your example python doesn't know it for me  
> **username_0**: when I print the variables of the nc file, it is there along y_range and z_range  

---

## 1 answers

### Answer by username_2 on 2022-06-15. Score: 5 | Accepted answer

It looks like you simply messed up the order of your min and max of the latitude. And you were looking at the wrong place (-44° instead of +44°).

```
# bathymetry
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

data = Dataset("ETOPO1_Bed_g_gmt4.grd",'r')

lon_range = data.variables['x_range'][:]
lat_range = data.variables['y_range'][:]
topo_range = data.variables['z_range'][:]
spacing = data.variables['spacing'][:]
dimension = data.variables['dimension'][:]
z = data.variables['z'][:]
lon_num = dimension[0]
lat_num = dimension[1]

lon = np.linspace(lon_range[0],lon_range[1],dimension[0])
lat = np.linspace(lat_range[0],lat_range[1],dimension[1])

topo = np.reshape(z, (lat_num, lon_num))

# Select region
lon_idx_min = np.abs(lon-10).argmin() # I increased the margin a little
lon_idx_max = np.abs(lon-20).argmin() # to better see where the region is
lat_idx_min = np.abs(lat+48).argmin()
lat_idx_max = np.abs(lat+40).argmin()

plt.imshow(topo[lat_idx_min:lat_idx_max, lon_idx_min:lon_idx_max], vmax=0)

```

By using np.linspace you can remove the loop from your code and speed it up significantly.
There is no need to create a meshgrid from lat and lon.

> **Comments:**
> **username_3**: greetings fellow meteorologist  

---

