remapbil: add support for healpix target grids

src/remap.h

@@ -218,7 +218,7 @@ std::pair<double, double> remap_find_weights(const LonLatPoint &llpoint, const d
 int rect_grid_search(size_t &ii, size_t &jj, double x, double y, size_t nxm, size_t nym, const Varray<double> &xm,
                      const Varray<double> &ym);
 
-LonLatPoint remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add);
+LonLatPoint remapgrid_get_lonlat(RemapGrid *grid, size_t index);
 
 void remap_check_area(size_t grid_size, const Varray<double> &cell_area, const char *name);
 

src/remaplib.cc

@@ -169,23 +169,28 @@ boundboxFromCenter(bool lonIsCyclic, size_t size, size_t nx, size_t ny, const Va
 }
 
 LonLatPoint
-remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add)
+remapgrid_get_lonlat(RemapGrid *grid, size_t index)
 {
   double lon, lat;
 
   if (grid->type == RemapGridType::Reg2D)
     {
       auto nx = grid->dims[0];
-      auto iy = cell_add / nx;
-      auto ix = cell_add - iy * nx;
+      auto iy = index / nx;
+      auto ix = index - iy * nx;
       lat = grid->reg2d_center_lat[iy];
       lon = grid->reg2d_center_lon[ix];
       if (lon < 0) lon += PI2;
     }
+  else if (grid->type == RemapGridType::HealPix)
+    {
+      hp_index_to_lonlat(grid->order, grid->nside, index, &lon, &lat);
+      //if (lon < 0) lon += PI2;
+    }
   else
     {
-      lat = grid->cell_center_lat[cell_add];
-      lon = grid->cell_center_lon[cell_add];
+      lat = grid->cell_center_lat[index];
+      lon = grid->cell_center_lon[index];
     }
 
   return LonLatPoint(lon, lat);
@@ -431,15 +436,16 @@ remap_define_grid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *
   int gridID_gme = -1;
 
   auto gridtype = gridInqType(grid.gridID);
+  auto isHealpixGrid = (grid.type == RemapGridType::HealPix);
 
-  if (grid.type == RemapGridType::HealPix)
+  if (isHealpixGrid)
     {
       auto [nside, order] = cdo_get_healpix_params(gridID);
       grid.nside = nside;
       grid.order = order;
     }
 
-  if (gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR && grid.type != RemapGridType::HealPix)
+  if (gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR && !isHealpixGrid)
     {
       constexpr auto withBounds = true;
       if (gridtype == GRID_GME)
@@ -464,9 +470,9 @@ remap_define_grid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *
 
   grid.size = gridInqSize(gridID);
 
-  grid.dims[0] = (grid.type == RemapGridType::HealPix) ? grid.size : gridInqXsize(gridID);
+  grid.dims[0] = isHealpixGrid ? grid.size : gridInqXsize(gridID);
   grid.dims[1] = gridInqYsize(gridID);
-  if (gridtype != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_UNSTRUCTURED && grid.type != RemapGridType::HealPix)
+  if (gridtype != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_UNSTRUCTURED && !isHealpixGrid)
     {
       if (grid.dims[0] == 0) cdo_abort("%s grid without longitude coordinates!", gridNamePtr(gridtype));
       if (grid.dims[1] == 0) cdo_abort("%s grid without latitude coordinates!", gridNamePtr(gridtype));
@@ -474,7 +480,7 @@ remap_define_grid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *
 
   grid.isCyclic = (gridIsCircular(gridID) > 0);
 
-  grid.rank = (gridInqType(gridID) == GRID_UNSTRUCTURED || grid.type == RemapGridType::HealPix) ? 1 : 2;
+  grid.rank = (gridInqType(gridID) == GRID_UNSTRUCTURED || isHealpixGrid) ? 1 : 2;
 
   grid.num_cell_corners = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? gridInqNvertex(gridID) : 4;
 
@@ -485,7 +491,7 @@ remap_define_grid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *
 
   if (!remap_write_remap && grid.type == RemapGridType::Reg2D) return;
 
-  if (grid.type != RemapGridType::HealPix && !gridHasCoordinates(gridID))
+  if (!isHealpixGrid && !gridHasCoordinates(gridID))
     cdo_abort("%s grid cell center coordinates missing!", txt);
 
   gridInqXvals(gridID, grid.cell_center_lon.data());
@@ -502,7 +508,7 @@ remap_define_grid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *
   if (gridInqType(grid.gridID) == GRID_GME) gridInqMaskGME(gridID_gme, &grid.vgpm[0]);
 
   // Convert lat/lon units if required
-  if (grid.type != RemapGridType::HealPix)
+  if (!isHealpixGrid)
     {
       cdo_grid_to_radian(gridID, CDI_XAXIS, grid.size, grid.cell_center_lon.data(), "grid center lon");
       cdo_grid_to_radian(gridID, CDI_YAXIS, grid.size, grid.cell_center_lat.data(), "grid center lat");
@@ -757,6 +763,12 @@ remap_init_grids(RemapMethod mapType, bool doExtrapolate, int gridID1, RemapGrid
         tgtGrid.type = RemapGridType::Reg2D;
     }
 
+  if (!remap_gen_weights && is_healpix_grid(gridID2))
+    {
+      if (mapType == RemapMethod::BILINEAR)
+        tgtGrid.type = RemapGridType::HealPix;
+    }
+
   srcGrid.doExtrapolate = doExtrapolate;
 
   if (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)