Converting from 2D lat/lon matrix into 1D lat/lon array

user_16173560 0
Viewed 303

I am working with a netcdf file, that has no coordinates. My lat/lon values are sotred in the variables in a form of a matrix, in a way lat (x,y) and lon(x,y). My goal here is to extract a lat and a lon 1D array to assign it as coordinates, as these must be a 1D array.

Here is how the dataset originally looks like:

<xarray.Dataset>
Dimensions:             (y: 10980, x: 10980)
Dimensions without coordinates: y, x
Data variables: (12/20)
    lon                 (y, x) float32 ...
    lat                 (y, x) float32 ...

For example, the lat variable looks like this:

<xarray.DataArray 'lat' (y: 10980, x: 10980)>
array([[41.52681 , 41.52681 , 41.526814, ..., 41.54671 , 41.54671 , 41.546715],
       [41.52672 , 41.526722, 41.526722, ..., 41.54662 , 41.546623, 41.546623],
       [41.52663 , 41.52663 , 41.526634, ..., 41.54653 , 41.54653 , 41.54653 ],
       ...,
       [40.538834, 40.538837, 40.538837, ..., 40.55806 , 40.55806 , 40.558064],
       [40.538746, 40.538746, 40.53875 , ..., 40.55797 , 40.557972, 40.557972],
       [40.538654, 40.53866 , 40.53866 , ..., 40.55788 , 40.55788 , 40.55788 ]],
      dtype=float32)
Dimensions without coordinates: y, x
Attributes:
    parameter:      lat
    standard_name:  latitude
    long_name:      latitude
    units:          degree_north

So in order to convert both variables into 1D array I do the following:

#First I open the dataset
file_to_input = 'landsat.nc'
nc1 = xr.open_dataset(file_to_input)

#Then I take the y axis from lat:
lati = nc1.lat[:,0]
#And the x axis from lon:
long = nc1.lon[0,:]

#To then assign them as 1D array to the dataset:
nc1 = nc1.assign_coords({'x':long,'y':lati})
nc1 = nc1.rio.set_spatial_dims('x', 'y')

#I set the proper CRS for the varaible to export (EPSG: 32631):
nc1var = nc1['ndci'].rio.set_crs("epsg:32631")

#And then I export it as a geotiff:
nc1var.rio.to_raster('ndci.tiff')

So far so good. The problem comes when I visualize the exported geotiff, the geotiff is SLIGHTLY shifted. In the image below you can appreciate a small shift down of the geotiff with respect to the basemap. I tried using this method with other good working tiffs and the same shift happens, so I assume it has to do witht he way I am changing from 2D lat lon to 1D array.

enter image description here

I think it can be achieved with pyproj transformer or smthing like that, but I have no idea how to use that with lat/lon 2D grid into 1D. Any help would be very appreciated!!

UPDATE

Download the dataset here

(.rar - 241MB, .nc file - 1.34GB)

python-xarray netcdf coordinates matrix python
1 Answers

I will propose following code to get data to regular grid using nearest neighbor interpolation. Other methods could be used for the interpolation, this is just the easiest and not changing original data much.

Also notice that I am reducing the original resolution for memory efficiency, this part should be otherwise removed from the script.

Here is the code:

#!/usr/bin/env ipython
# ---------------------
import numpy as np
from scipy.interpolate import griddata
from pylab import pcolormesh, show
from netCDF4 import Dataset
# --------------------------------------------
def nc_varget(fin,vin):
    with Dataset(fin) as f: return f.variables[vin][:];
# --------------------------------------------
fin = 'S2A_reduced_dataset.nc' 
xin = nc_varget(fin,'lon');
yin = nc_varget(fin,'lat');
zin = nc_varget(fin,'ndci');
# --------------------------------------------
# I will reduce the size further to test the code (matrices with 10k X 10k points are too big for my laptop and testing purposes)
xskip, yskip = 25,25
xin = xin[::yskip,::xskip]
yin = yin[::yskip,::xskip]
zin = zin[::yskip,::xskip]
# --------------------------
# let us take some info from original coordinates:
x0,x1,dx = np.min(xin),np.max(xin),np.abs(np.mean(np.diff(xin)))
y0,y1,dy = np.min(yin),np.max(yin),np.abs(np.mean(np.diff(yin.T)))
# --------------------------
# let us make new (regular) coordinates:
xout = np.arange(x0,x1+dx,dx)
yout = np.arange(y0,y1+dx,dy)
# --------------------------
xm,ym = np.meshgrid(xout,yout)
zo = griddata((xin.flatten(),yin.flatten()),zin.flatten(),(xm,ym),'nearest')
# ---------------------------------------------------------------------------
# let us save results as netCDF:
import xarray as xr
df = xr.DataArray(zo,dims=['lat','lon'])
df.coords['lon'] = xout
df.coords['lat'] = yout
df.to_netcdf('test.nc')
user_19133670
Related