4
$\begingroup$

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...)

$\endgroup$
2
  • $\begingroup$ hi, where does x_range come from, because when i try to rebuild your example python doesn't know it for me $\endgroup$
    – Weiss
    Commented Jun 22, 2022 at 10:46
  • $\begingroup$ when I print the variables of the nc file, it is there along y_range and z_range $\endgroup$
    – shamalaia
    Commented Jun 22, 2022 at 23:28

1 Answer 1

5
+50
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ greetings fellow meteorologist $\endgroup$ Commented Jun 15, 2022 at 10:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.