Need help interpolating data into specific points in netcdf files

Viewed 40

I need help interpolating data onto nodes in a series of netcdf files which have oceanographic data. The files have 3 dimensions (latitude (12), longitude (20), time (24)) and variables (u, v) . Of all the data points, there are four nodes which do not have current velocity data (u, v) although they should have data, as they are at sea but register as land. I am trying to interpolate data onto these nodes but have no idea how.

EDIT: The four points with missing data are already in the coordinates but are have NaN values. The other points would keep the original data.

I am OK with Pandas but I know that this probably requires numpy and/or xarray and I am not well versed in either. I cannot get to nodes using the coordinates to interpolate the data I require. Can this be done at all?

print(data)
<xarray.Dataset>
Dimensions:    (latitude: 12, time: 24, longitude: 20)
Coordinates:
* latitude   (latitude) float32 40.92 41.0 41.08 41.17 ... 
41.67 41.75 41.83
* time       (time) datetime64[ns] 2017-03-03T00:30:00 ... 
2017-03-03T23:30:00
* longitude  (longitude) float32 1.417 1.5 1.583 1.667 ... 
2.833 2.917 3.0
Data variables:
    vo         (time, latitude, longitude) float32 ...
    uo         (time, latitude, longitude) float32 ...
    thetao     (time, latitude, longitude) float32 ...
    zos        (time, latitude, longitude) float32 ...
Attributes: (12/22)
Conventions:                        CF-1.0
source:                             CMEMS IBI-MFC...


print(data.latitude)
<xarray.DataArray 'latitude' (latitude: 12)>
array([40.916668, 41.      , 41.083332, 41.166668, 41.25    , 
41.333332, 41.416668, 41.5     , 41.583332, 41.666668, 41.75    , 41.833332], dtype=float32)
Coordinates:
* latitude  (latitude) float32 40.92 41.0 41.08 41.17 ... 41.67 
41.75 41.83
Attributes:
standard_name:        latitude
long_name:            Latitude
units:                degrees_north
axis:                 Y
unit_long:            Degrees North
step:                 0.08333f
_CoordinateAxisType:  Lat
_ChunkSizes:          361
valid_min:            40.916668
valid_max:            41.833332

print(data.longitude)
<xarray.DataArray 'longitude' (longitude: 20)>
array([1.416666, 1.499999, 1.583333, 1.666666, 1.749999, 1.833333, 1.916666, 1.999999, 2.083333, 2.166666, 2.249999, 2.333333, 2.416666, 2.499999,2.583333, 2.666666, 2.749999, 2.833333, 2.916666, 2.999999],
  dtype=float32)
Coordinates:
* longitude  (longitude) float32 1.417 1.5 1.583 1.667 ... 2.833 2.917 3.0
Attributes:
standard_name:        longitude
long_name:            Longitude
units:                degrees_east
axis:                 X
unit_long:            Degrees East
step:                 0.08333f
_CoordinateAxisType:  Lon
_ChunkSizes:          289
valid_min:            1.4166658
valid_max:            2.999999
1 Answers

The goal of the question is to in-fill cells, which are on land/missing in the raw files, but should really be in the sea. In some cases, you might want something sophisticated to do this. For example, if there was a sharp coastal gradient.

But the easiest way to solve it is to use nearest neighbour to replace missing values with the nearest neighbour. That will of course replace more than you need. So, you will then need to apply some kind of land-sea mask to your data. The workflow below, using my package nctoolkit, should do the job

import nctoolkit as nc
# read in the file and set missing values to nn
ds = nc.open_data("infile.nc")
ds.fill_na()
# create the mask to apply. This should only have one time step
# I'm going to assume in this case that it is a file with temperature that has the correct land-sea division
ds_mask = nc.open_data("mask.nc")
# Ensure sea values are 1. Land values should be nan
ds_mask.compare(">-1000")
# multiply the dataset by the mask to set land values to missing
ds.multiply(ds_mask)
# plot the results
ds.plot()
Related