I have a raster dataset (NetCDF) but I can't find the correct projection for it. I'm reaching out so maybe someone can give me some indication of the right projection. Currently the project CRS is in ETRS89 (EPSG:4258), as shown in the image attached. I tried many CRS and both raster and vector layers just don't line up: there always seems to be N-S shift. The data is sourced from this link: https://thredds.socib.es/lw4nc2/index.html?m=wmop_surface Any insight?
The problem is that the data is not defined on a regularly spaced grid.
In a NetCDF file the coordinates of a grid do not have to be evenly spaced, they are set by the "dimension" in lat and long. A 3x3 grid could have latitude dimension of [1,2,3] which gives regular-spaced 1 degree cells, or [1, 2.2, 3] which would be one row of 1.2 degree height and one row of 0.8 degree height.
If your software can't handle irregular grids (QGIS cant, which is what I'm using) then it takes the min and max and assumes regular grid latitudes. In the simple case above that would mean a true latitude of [1,2.2,3] would be interpreted as [1,2,3]. Any grids that should be at latitude=2.2 will be drawn at latitude 2 instead.
This looks like what is happening here, and QGIS warns me.
Warning 1: Latitude grid not spaced evenly.
Setting projection for grid spacing is within 0.1 degrees threshold.
The most northerly point (around Genoa) is correct, as is the most southerly (E of Melila). Between those latitudes there's a shift. Its not as simple as moving the whole raster up.
I can inspect the grid coordinate values in R using its ncdf4 package. For example, the longitude values are all 0.02380952 degrees apart, so its a regular longitude grid:
> head(lon$vals) # get first few longitudes
[1] -5.800000 -5.776190 -5.752381 -5.728571 -5.704762 -5.680952
Whats the smallest and biggest difference between successive longitudes?
> range(diff(lon$vals))
[1] 0.02380952 0.02380952
But for latitude:
> range(diff(lat$vals))
[1] 0.01692497 0.01952743
So the smallest latitude bands are 0.0169 degrees and the largest are 0.0195 degrees. This, coupled with the GIS fitting the grid between the extents, gives the shape you see.
The fix is to use a system that can cope with "curvilinear" grids, which I think the only example I know is the stars package in R (but not sure) and some other specialist remote sensing software.
The other solution is to treat the data as (x, y, z) triples and re-interpolate to a regular grid in whatever coordinate system you want - this doesn't even have to be lat-long you could use some cartesian grid (with care).
