🔀 Merge branch 'egu2019_poster_figs' (#136)

Closes #136 Reproduce figures in EGU2019 DeepBedMap poster.

Pipfile

@@ -26,7 +26,6 @@ pandas = "==0.24.2"
 pyproj = "==2.1.0"
 quilt = "==2.9.14"
 rasterio = "==1.0.13"
-requests = "==2.21.0"
 scikit-image = "==0.14.2"
 shapely = "==1.7a1"
 toolz = "==0.9.0"

data_list.yml

@@ -173,3 +173,16 @@
        filename: WISE_ISODYN_RadarByFlight_ASCII.zip
        url: https://ramadda.data.bas.ac.uk/repository/entry/get/WISE_ISODYN_RadarByFlight_ASCII.zip?entryid=synth%3A59e5a6f5-e67d-4a05-99af-30f656569401%3AL1dJU0VfSVNPRFlOX1JhZGFyQnlGbGlnaHRfQVNDSUkuemlw
        sha256: dfb20a4ff133f361f64b5730e7593a46cb3d599ae62b0dc874685794990c8d91
+-
+  citekey: Fahnestock2019LISA
+  folder: misc
+  location: Antarctica
+  resolution: 750m
+  doi:
+    dataset: "https://doi.org/10.7265/nxpc-e997"
+    literature: "https://doi.org/10.1016/j.rse.2015.11.023"
+  files:
+    -
+      filename: lisa750_2013182_2017120_0000_0400_vv_v1.tif
+      url: ftp://ftp.nsidc.org/pub/DATASETS/nsidc0733_landsat_ice_speed_v01/LISA750/lisa750_2013182_2017120_0000_0400_v1.tgz
+      sha256: 99eb934702305e2d27afa20dfc211c1d45ed762727bda29e0f251976a1877a92

data_prep.py

@@ -1,12 +1,12 @@
 # ---
 # jupyter:
 #   jupytext:
-#     formats: ipynb,py:percent
+#     formats: ipynb,py:hydrogen
 #     text_representation:
 #       extension: .py
-#       format_name: percent
+#       format_name: hydrogen
 #       format_version: '1.2'
-#       jupytext_version: 1.0.1
+#       jupytext_version: 1.0.3
 #   kernelspec:
 #     display_name: deepbedmap
 #     language: python
@@ -29,8 +29,9 @@
 import os
 import shutil
 import sys
+import tarfile
+import urllib
 
-import requests
 import tqdm
 import yaml
 
@@ -63,31 +64,53 @@
 # %%
 def download_to_path(path: str, url: str):
     r"""
-    Download from a url to a path
+    Download from a HTTP or FTP url to a filepath.
 
-    >>> download_to_path(path="highres/Data_20171204_02.csv",
-    ...                  url="https://data.cresis.ku.edu/data/rds/2017_Antarctica_Basler/csv_good/Data_20171204_02.csv")
-    <Response [200]>
+    >>> d = download_to_path(
+    ...    path="highres/Data_20171204_02.csv",
+    ...    url="ftp://data.cresis.ku.edu/data/rds/2017_Antarctica_Basler/csv_good/Data_20171204_02.csv",
+    ... )
     >>> open("highres/Data_20171204_02.csv").readlines()
     ['LAT,LON,UTCTIMESOD,THICK,ELEVATION,FRAME,SURFACE,BOTTOM,QUALITY\n']
     >>> os.remove(path="highres/Data_20171204_02.csv")
     """
-    # if not os.path.exists(path=path):
-    r = requests.get(url=url, stream=True)
-    with open(file=path, mode="wb") as fd:
-        for chunk in r.iter_content(chunk_size=1024):
-            fd.write(chunk)
-    return r
+
+    folder, filename = os.path.split(p=path)
+    downloaded_filename = os.path.basename(urllib.parse.urlparse(url=url).path)
+
+    # Download file using URL first
+    if not os.path.exists(os.path.join(folder, downloaded_filename)):
+        r = urllib.request.urlretrieve(
+            url=url, filename=os.path.join(folder, downloaded_filename)
+        )
+
+    # If downloaded file is not the final file (e.g. file is in an archive),
+    # then extract the file from the archive!
+    if filename != downloaded_filename:
+        # Extract tar.gz archive file
+        if downloaded_filename.endswith(("tgz", "tar.gz")):
+            try:
+                archive = tarfile.open(name=f"{folder}/{downloaded_filename}")
+                archive.extract(member=filename, path=folder)
+            except:
+                raise
+        else:
+            raise ValueError(
+                f"Unsupported archive format for downloaded file: {downloaded_filename}"
+            )
+
+    return os.path.exists(path=path)
 
 
 # %%
 def check_sha256(path: str):
     """
     Returns SHA256 checksum of a file
 
-    >>> download_to_path(path="highres/Data_20171204_02.csv",
-    ...                  url="https://data.cresis.ku.edu/data/rds/2017_Antarctica_Basler/csv_good/Data_20171204_02.csv")
-    <Response [200]>
+    >>> d = download_to_path(
+    ...    path="highres/Data_20171204_02.csv",
+    ...    url="https://data.cresis.ku.edu/data/rds/2017_Antarctica_Basler/csv_good/Data_20171204_02.csv",
+    ... )
     >>> check_sha256("highres/Data_20171204_02.csv")
     '53cef7a0d28ff92b30367514f27e888efbc32b1bda929981b371d2e00d4c671b'
     >>> os.remove(path="highres/Data_20171204_02.csv")
@@ -118,7 +141,7 @@ def parse_datalist(
     assert yaml_file.endswith((".yml", ".yaml"))
 
     with open(file=yaml_file, mode="r") as yml:
-        y = yaml.load(stream=yml)
+        y = yaml.safe_load(stream=yml)
 
     datalist = pd.io.json.json_normalize(
         data=y, record_path=record_path, meta=schema, sep="_"
@@ -192,7 +215,7 @@ def parse_datalist(
     rasterio.plot.show(source=raster_source, cmap="BrBG_r")
 
 # %% [markdown]
-# ### Download miscellaneous data (e.g. [REMA](https://doi.org/10.7910/DVN/SAIK8B), [MEaSUREs Ice Flow](https://doi.org/10.5067/OC7B04ZM9G6Q))
+# ### Download miscellaneous data (e.g. [REMA](https://doi.org/10.7910/DVN/SAIK8B), [MEaSUREs Ice Flow](https://doi.org/10.5067/D7GK8F5J8M8R), [LISA](https://doi.org/10.7265/nxpc-e997))
 
 # %%
 for dataset in dataframe.query(expr="folder == 'misc'").itertuples():
@@ -234,9 +257,10 @@ def ascii_to_xyz(pipeline_file: str) -> pd.DataFrame:
     a JSON Pipeline file similar to the one used by PDAL.
 
     >>> os.makedirs(name="/tmp/highres", exist_ok=True)
-    >>> download_to_path(path="/tmp/highres/2011_Antarctica_TO.csv",
-    ...                  url="https://data.cresis.ku.edu/data/rds/2011_Antarctica_TO/csv_good/2011_Antarctica_TO.csv")
-    <Response [200]>
+    >>> d = download_to_path(
+    ...    path="/tmp/highres/2011_Antarctica_TO.csv",
+    ...    url="https://data.cresis.ku.edu/data/rds/2011_Antarctica_TO/csv_good/2011_Antarctica_TO.csv",
+    ... )
     >>> _ = shutil.copy(src="highres/20xx_Antarctica_TO.json", dst="/tmp/highres")
     >>> df = ascii_to_xyz(pipeline_file="/tmp/highres/20xx_Antarctica_TO.json")
     >>> df.head(2)
@@ -538,7 +562,7 @@ def selective_tile(
     ...    filepath="/tmp/tmp_st.nc",
     ...    window_bounds=[(1.0, 4.0, 3.0, 6.0), (2.0, 5.0, 4.0, 7.0)],
     ... )
-    Tiling: /tmp/tmp_st.nc
+    Tiling: /tmp/tmp_st.nc ... done!
     array([[[[0.18485446, 0.96958464],
              [0.4951769 , 0.03438852]]],
     <BLANKLINE>
@@ -550,7 +574,7 @@ def selective_tile(
     array_list = []
 
     with rasterio.open(filepath) as dataset:
-        print(f"Tiling: {filepath}")
+        print(f"Tiling: {filepath} ... ", end="")
         for window_bound in window_bounds:
 
             if padding > 0:
@@ -595,15 +619,18 @@ def selective_tile(
                     np.copyto(
                         dst=array, src=array2, where=array.mask
                     )  # fill in gaps where mask is True
+
+                    # assert not array.mask.any()  # ensure no NAN values after gapfill
                 else:
                     plt.imshow(array.data[0, :, :])
                     plt.show()
-                    raise ValueError(
-                        f"Tile has missing data, try passing in gapfill_raster_filepath"
+                    print(
+                        f"WARN: Tile has missing data, try passing in gapfill_raster_filepath"
                     )
 
             # assert array.shape[0] == array.shape[1]  # check that height==width
             array_list.append(array.data.astype(dtype=np.float32))
+        print("done!")
 
     return np.stack(arrays=array_list)
 
@@ -649,12 +676,21 @@ def selective_tile(
 )
 print(rema.shape, rema.dtype)
 
+# %%
+## Custom processing for LISA to standardize units with MEASURES Ice Velocity
+# Convert units from metres/day to metres/year by multiplying 1st band by 365.25
+!rio calc "(* 365.25 (read 1))" misc/lisa750_2013182_2017120_0000_0400_vv_v1.tif misc/lisa750_2013182_2017120_0000_0400_vv_v1_myr.tif
+# Set NODATA mask where pixels are 36159.75 = 99 * 365.25
+!rio edit-info misc/lisa750_2013182_2017120_0000_0400_vv_v1_myr.tif --nodata 36159.75
+!rio info misc/lisa750_2013182_2017120_0000_0400_vv_v1_myr.tif
+
 # %%
 measuresiceflow = selective_tile(
     filepath="misc/MEaSUREs_IceFlowSpeed_450m.tif",
     window_bounds=window_bounds_concat,
     padding=1000,
     out_shape=(20, 20),
+    # gapfill_raster_filepath="misc/lisa750_2013182_2017120_0000_0400_vv_v1_myr.tif",
 )
 print(measuresiceflow.shape, measuresiceflow.dtype)