Selecting the Indian Ocean region (which includes the seam) from ACCESS-OM ocean model output

COSIMA Ocean
1

I would like to select Indian Ocean 40E-100E, 20S-20N

idir = '/g/data/p66/mrd599/cq880/history/ocn/'
fname = idir + 'ocean-2d-surface_pot_temp-1-monthly-mean-ym_03*.nc'

ds = xr.open_mfdataset(fname, combine='by_coords')
ds.surface_pot_temp
lat_range = (-20, 20)
lon_range = (40, 100)
ds2 = ds.assign_coords({"xt_ocean": (ds["xt_ocean"] + 360)  % 360}) # <- previous suggested but not working

lon_indices = np.where((ds2.xt_ocean >= lon_range[0]) & (ds2.xt_ocean <= lon_range[1]))[0]
lat_indices = np.where((ds.yt_ocean >= lat_range[0]) & (ds.yt_ocean <= lat_range[1]))[0]
print(lon_indices)

The output does not make sense:

[   0    1    2    3    4    5    6    7    8    9   10   11   12   13
   14   15   16   17   18   19   20   21   22   23   24   25   26   27
   28   29   30   31   32   33   34   35   36   37   38   39   40   41
   42   43   44   45   46   47   48   49   50   51   52   53   54   55
   56   57   58   59   60   61   62   63   64   65   66   67   68   69
   70   71   72   73   74   75   76   77   78   79 1280 1281 1282 1283
 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297
 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311
 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325
 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339
 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353
 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367
 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381
 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395
 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409
 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423
 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437
 1438 1439]
2

This just resets the coordinate values to be in the range 0-360 deg, but the Indian Ocean is still split across the edge of the array. There are a number of ways of dealing with slicing periodic coordinates across array edges, but the simplest here is to sort by xt_ocean prior to slicing out the Indian Ocean. E.g., picking up from your code:

ds_indian_ocean = ds2.sortby("xt_ocean").sel(
    xt_ocean=slice(40, 100),
    yt_ocean=slice(-20, 20),
)
2 Likes
3

Many thanks!! I will give a try now.