Merge pull request #215 from nsidc/fix-warp-hack

Fix warp hack

CHANGELOG.md

@@ -1,3 +1,8 @@
+# v1.0.0rc2 (2021-02-11)
+
+- Fixed mistake in raster processing code that resulted in GDAL unnecessarily
+  requesting ~22GB of diskspace for reprojection.
+
 # v1.0.0rc1 (2021-02-10)
 
 - Add styles for World Magnetic Model layers.

qgreenland/__init__.py

@@ -1 +1 @@
-__version__ = 'v1.0.0rc1'
+__version__ = 'v1.0.0rc2'

qgreenland/config/layers.yml

@@ -1227,6 +1227,8 @@
   decompress_kwargs:
     extract_files:
       - 'magnitudes.nc'
+  warp_kwargs:
+    ignore_output_bounds_hack: True
 
 - &racmo_1km
   id: racmo_precip
@@ -1252,6 +1254,8 @@
       - 'precip.1958-2019.BN_RACMO2.3p2_FGRN055_1km.YY-mean.nc'
   unzip_kwargs:
     input_filename: 'RACMO_QGreenland_Jan2021.zip'
+  warp_kwargs:
+    ignore_output_bounds_hack: True
 
 - <<: *racmo_1km
   id: racmo_snowfall
@@ -1357,6 +1361,8 @@
   decompress_kwargs:
     extract_files:
       - 'Icemask_Topo_Iceclasses_lon_lat_average_1km_GrIS.nc'
+  warp_kwargs:
+    ignore_output_bounds_hack: True
 
 - <<: *racmo_1km_masks
   id: racmo_grounded_ice
@@ -1390,6 +1396,7 @@
       - 'Icemask_Topo_Iceclasses_lon_lat_average_1km_GrIS.nc'
   warp_kwargs:
     srcNodata: 9999.0
+    ignore_output_bounds_hack: True
 
 ##############
 # Glaciology #
@@ -1635,9 +1642,7 @@
     extract_files:
       - 'ICESat1_ICESat2_mass_change/gris.tif'
   warp_kwargs:
-    # This dataset is in 3413 although the metadata does not state this
-    # explicitly.
-    srcSRS: 'EPSG:3413'
+    ignore_output_bounds_hack: True
 
 - &cci_surface_elevation_change
   id: surface_elevation_change_1992_1996

qgreenland/util/raster.py

@@ -67,14 +67,17 @@ def warp_raster(inp_path, out_path, *, layer_cfg, warp_kwargs=None):
                            'No projection automatically detected and '
                            'none explicitly provided.')
 
-    # gdal.Warp(out_path, inp_path, **warp_kwargs)
-    _gdalwarp_cut_hack(
-        out_path, inp_path,
-        layer_cfg=layer_cfg, warp_kwargs=warp_kwargs, src_srs_str=srs_str
-    )
+    ignore_output_bounds_hack = warp_kwargs.pop('ignore_output_bounds_hack', False)
+    if ignore_output_bounds_hack:
+        _gdalwarp(out_path, inp_path, **warp_kwargs)
+    else:
+        _gdalwarp_cut_hack(
+            out_path, inp_path,
+            layer_cfg=layer_cfg, warp_kwargs=warp_kwargs
+        )
 
 
-def _gdalwarp_cut_hack(out_path, inp_path, *, layer_cfg, warp_kwargs, src_srs_str):
+def _gdalwarp_cut_hack(out_path, inp_path, *, layer_cfg, warp_kwargs):
     """Hack for an issue with some gdal cutlines.
 
     <mailing list link here>
@@ -101,22 +104,10 @@ def _gdalwarp_cut_hack(out_path, inp_path, *, layer_cfg, warp_kwargs, src_srs_st
     # don't, e.g. compress twice.
     step2_keys = ['cutlineDSName', 'cropToCutline', 'creationOptions']
 
-    ignore_output_bounds_hack = warp_kwargs.pop('ignore_output_bounds_hack', False)
-
     # Step 1 needs to subset for this to work (outputBounds == `-te`).
     step1_kwargs = {k: v for k, v in warp_kwargs.items() if k not in step2_keys}
 
-    # TODO: this hack is not perfect!! The `src_srs_str` may report WKT params
-    # equivilent to EPSG:3413 without it actually being the string `EPSG:3413`,
-    # which is what our project CRS is defined as.
-    # HACK ANOTHER TERRIBLE HACK: reprojection and clipping should really be two
-    # separate steps!!!
-    # if the source dataset is already in the project's CRS, do not use
-    # `outputBounds` to limit the raster's extent. If the boundary's bbox does
-    # not exactly align with the source grid, the resulting grid will be
-    # spatially offset from the source.
-    if src_srs_str != CONFIG['project']['crs'] and not ignore_output_bounds_hack:
-        step1_kwargs['outputBounds'] = layer_cfg['boundary']['bbox']
+    step1_kwargs['outputBounds'] = layer_cfg['boundary']['bbox']
 
     # Step 2 actually does the shape-based cut as a separate step, to avoid
     # errors.