4

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

shamalaia
  • 529
  • 2
  • 14

1 Answers1

5

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.

smichel
  • 380
  • 1
  • 4