I’m new here and would like to thank anyone that takes the time to respond !
I am trying to compare two netcdf’s from different sources. The first one (A) is on a longitude/latitude grid.
The second one (B) has a “center” longitude and latitude, but then the coordinates are x and y, the horizontal and vertical distances in km to that center point.
I would like to map the bias between the two. So I’m trying to convert B into lon/lat coordinates, then interpolate to A’s grid.
I also have the longitude and latitude of each grid point of B. So lon_B and lat_B are two-dimensional arrays too.
As a first step : how do I make B_new, the new data array with lon/lat coordinates ? (there would be some nan values since the lon/lat grid wouldn’t coincide exactly with the x/y grid, but I don’t see another way to move forward)
I have tried looking into pyproj tools, but I think I’m missing something in my comprehension of how I can create a new array with the old array’s data and a new coordinate system that doesn’t coincide with the old one.
Edit : following @Patrick’s recommendation, here are the ncdump headers for A and B :
netcdf tg_ens_mean_0.1deg_reg_v29.0e {
dimensions:
latitude = 465 ;
longitude = 705 ;
time = UNLIMITED ; // (27028 currently)
variables:
double latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "Latitude values" ;
latitude:axis = "Y" ;
latitude:standard_name = "latitude" ;
double longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "Longitude values" ;
longitude:axis = "X" ;
longitude:standard_name = "longitude" ;
short tg(time, latitude, longitude) ;
tg:units = "Celsius" ;
tg:_FillValue = -9999s ;
tg:long_name = "mean temperature" ;
tg:scale_factor = 0.01f ;
tg:add_offset = 0.f ;
tg:standard_name = "air_temperature" ;
tg:cell_methods = "time: mean" ;
int time(time) ;
time:units = "days since 1950-01-01 00:00" ;
time:long_name = "Time in days" ;
time:calendar = "standard" ;
time:standard_name = "time" ;
time:cell_methods = "time: mean" ;
// global attributes:
:E-OBS_version = "29.0e" ;
:Conventions = "CF-1.4" ;
:References = "http://surfobs.climate.copernicus.eu/ dataaccess/access_eobs.php" ;
:history = "Fri Mar 22 09:55:59 2024: ncks --no-abc -d time,0,27027 /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg//tg_ensmean_master_rectime.nc /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg//tg_ens_mean_0.1deg_reg_v29.0e.ncnFri Mar 22 09:50:30 2024: ncks --no-abc --mk_rec_dmn time /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg//tg_ensmean_master.nc /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg//tg_ensmean_master_rectime.ncnFri Mar 22 09:32:14 2024: ncks -O -d time,0,27027 /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg/tg_ensmean_master_untilJan2024.nc /nobackup_1/users/besselaa/Data/Gridding/EOBSv29.0e/Grid_0.1deg/tg/tg_ensmean_master.nc" ;
:NCO = "netCDF Operators version 5.1.4 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
}
netcdf ICE.ERA5.EUc.TTz.2015 {
dimensions:
time = UNLIMITED ; // (365 currently)
ztqlev = 3 ;
y = 126 ;
x = 201 ;
variables:
float TTz(time, ztqlev, y, x) ;
TTz:units = "C" ;
TTz:long_name = "Temperature" ;
TTz:standard_name = "Temperature" ;
TTz:actual_range = -20.69143f, 20.23214f ;
TTz:cell_methods = "time: mean" ;
float time(time) ;
time:units = "DAYS since 2014-09-01 00:00:00" ;
time:long_name = "time" ;
time:standard_name = "time" ;
time:cell_methods = "time: mean" ;
float x(x) ;
x:units = "km" ;
x:long_name = "x" ;
x:standard_name = "x" ;
float y(y) ;
y:units = "km" ;
y:long_name = "y" ;
y:standard_name = "y" ;
float ztqlev(ztqlev) ;
ztqlev:units = "m" ;
ztqlev:long_name = "ztqlev" ;
ztqlev:standard_name = "ztqlev" ;
// global attributes:
:title = "ICE - Exp: c01 - 20150101" ;
:institution = "ULg (Xavier Fettweis)" ;
:history = "Wed Jul 31 15:04:58 2024: ncks -v TTz /bettik/castelli/MARout/EUc/c01/2015/ICE.c01.2015.01.01-31.nc /bettik/castelli/data/MAR3.14/MAR-ERA5/EUc/ICE.ERA5.EUc.TTz.01.2015.ncnlibUN (2005.04.08) - Tue Jun 11 12:03:35 2024" ;
:netcdf = "4.8.0 of May 2 2023 12:43:48 $" ;
:NCO = "netCDF Operators version 5.0.6 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
}
4
Following your comment, the first step would be to project B’s grid in A’s CRS :
from pyproj import CRS, Transformer
import numpy as np
crs_b = CRS.from_proj4("+proj=stere +lat_0=44.6 +lon_0=7.2 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
x, y = ds_b['x'].values, ds_b['y'].values
xx, yy = np.meshgrid(x, y)
lon_b, lat_b = Transformer.from_crs(crs_from=crs_b, crs_to=4326, always_xy=True).transform(xx, yy)
Then, you could use interpolation routines ou regridding algorithms to project B’s data onto A’s grid, and then perform difference between both data sources.
You could also have a look to the rasterio package and the reproject method, it may suits your needs.
1