Support for dealing with rotated pole coordinates #3593
Description
Hello all,
First of all I'm going to introduce the necessity.
I'm working with NetCDF files whose latitude and longitude dimensions are rotated instead of, say, "normal". These files contain CORDEX climate data from global and regional climate models. The file metadata is as follows:
double rlon(rlon=210);
:standard_name = "grid_longitude";
:long_name = "longitude in rotated pole grid";
:units = "degrees";
:axis = "X";
:_ChunkSizes = 210U; // uint
double rlat(rlat=139);
:standard_name = "grid_latitude";
:long_name = "latitude in rotated pole grid";
:units = "degrees";
:axis = "Y";
:_ChunkSizes = 139U; // uint
char rotated_pole;
:grid_mapping_name = "rotated_latitude_longitude";
:grid_north_pole_latitude = 75.74; // double
:grid_north_pole_longitude = -246.01999999999998; // double
float sfcWindmax(time=1825, rlat=139, rlon=210);
:grid_mapping = "rotated_pole";
:_FillValue = 1.0E20f; // float
:missing_value = 1.0E20f; // float
:standard_name = "wind_speed";
:long_name = "Daily Maximum Near-Surface Wind Speed";
:units = "m s-1";
:coordinates = "lon lat height";
:cell_methods = "time: maximum";
:_ChunkSizes = 1U, 139U, 210U; // uint
double lon(rlat=139, rlon=210);
:standard_name = "longitude";
:long_name = "longitude";
:units = "degrees_east";
:_ChunkSizes = 139U, 210U; // uint
double lat(rlat=139, rlon=210);
:standard_name = "latitude";
:long_name = "latitude";
:units = "degrees_north";
:_ChunkSizes = 139U, 210U; // uint
The data variable of interest is "sfcWindmax" which is dimensioned along rotated coordinates "rlat" and "rlon". When reading it with netcdf-java library, I can extract all the data into an Array and therefore create a Geotrellis 210x139 Tile containing the data.
The problem arises when I try to interpret the extent of that tile, which, a priori, would be (, , , ), but those coordinates are rotated. The real extent of this NetCDF file is CAM region, (-127.8918101637883922,-19.9447726068228945,-19.0821144147941908,46.6000000000000085), whereas the rotated extent is (-52.800000000000004, -28.6, 39.15999999999999, 32.12). So, in order to treat this raster with real world EPSG:4326 coordinates, a coordinates transformation needs to be performed, in other words, a raster reprojection.
When showing a PNG of the tile with the rotated extent I get a perfectly rectangular tile:
By using QGIS (GDAL), I am able to reproject it to get the following raster:
Which is the raster I want to use in my process, since it has the extent with real coordinates, so that's the target raster I want to get with Geotrellis.
As a first attempt I have manually reprojected the raster. The data providers already have done that transformation task (from rotated to non rotated coordinates) for the user and give you another two-dimension variables "lon" and "lat" that contain the transformation of every pair (rlat, rlon). The target raster shown is a bit larger than the 210x139 original, it is approx 243x148. So one method could be creating an empty 243x148 tile and fill it by mapping every pair of (lat, lon) -> sfcWindmax real coordinates to its cell, with an extent (-127.8918101637883922,-19.9447726068228945,-19.0821144147941908,46.6000000000000085).
The result is somehow similar to the target but there are lots of cells with no value, since there are cells that are not mapped from any pair of coordinates:
So I have created a rudimentary algorithm to fill the blank cells by KNN or Bilinear sampling.
At this point, I already know that there is the method "reproject" on a tile, with some variations, for this need. The problem is that I cannot manage to get a result; every attempt of reprojecti via that method results in a tile filled with NaN.
The best result I can get is the following raster, which indeed it is CAM region but raster values and the extent do not match the expected:
So, please, has anybody dealt with this problem before?
I cannot use GDAL because this process is expected to perform in a EMR cluster.
I have already tried reprojecting via proj4j library but it lacks the Rotation Projection, which is the one I need to use.
Thank you.