diff --git a/antarctica_today/__init__.py b/antarctica_today/__init__.py
index e869860..624a501 100644
--- a/antarctica_today/__init__.py
+++ b/antarctica_today/__init__.py
@@ -1,4 +1,5 @@
 import os
+import sys
 
 # IMPORTANT: If we don't specify this setting, then the projection we want to use will
 # be replaced with another (and this warning will be printed)!
@@ -10,3 +11,11 @@
 #   used instead. To use the original CRS, set the OSR_USE_NON_DEPRECATED configuration
 #   option to NO.
 os.environ["OSR_USE_NON_DEPRECATED"] = "NO"
+
+
+# Ignore warnings by default, while still allowing users to change the behavior, e.g. by
+# upgrading them to exceptions.
+if not sys.warnoptions:
+    import warnings
+
+    warnings.simplefilter("ignore")
diff --git a/antarctica_today/generate_antarctica_today_map.py b/antarctica_today/generate_antarctica_today_map.py
index 877059e..dceede0 100644
--- a/antarctica_today/generate_antarctica_today_map.py
+++ b/antarctica_today/generate_antarctica_today_map.py
@@ -54,7 +54,8 @@ def main():
     """Do stuff I want to do here."""
 
     m = AT_map_generator(
-        fill_pole_hole=False, filter_out_error_swaths=True, verbose=True
+        fill_pole_hole=False,
+        filter_out_error_swaths=True,
     )
 
     for region in [
@@ -208,13 +209,6 @@ def read_and_parse_args():
         default=False,
         help="Omit the legend. Default if not set: include a legend.",
     )
-    parser.add_argument(
-        "--verbose",
-        "-v",
-        action="store_true",
-        default=False,
-        help="Increase output verbosity.",
-    )
 
     return parser.parse_args()
 
@@ -231,7 +225,6 @@ def __init__(
         melt_array_picklefile=model_results_picklefile,
         fill_pole_hole=True,
         filter_out_error_swaths=True,
-        verbose=True,
     ):
         """Initialize the class."""
         self.melt_array_picklefile = melt_array_picklefile
@@ -893,7 +886,7 @@ def _strip_empty_image_border(self, filename):
             return
         # svgclip.py isn't working... can't seem to resolve the Rsvg namespace.
         # svgclip.clip(filename, filename, margin=0)
-        # logger.debug("Trimmed {filename}.")
+        # logger.debug(f"Trimmed {filename}")
 
         else:
             bg = PIL.Image.new(im.mode, im.size, im.getpixel((0, 0)))
@@ -903,7 +896,7 @@ def _strip_empty_image_border(self, filename):
             if bbox:
                 im2 = im.crop(bbox)
                 im2.save(filename)
-                logger.debug("Trimmed {filename}.")
+                logger.debug(f"Trimmed {filename}")
 
         return
 
@@ -1790,7 +1783,6 @@ def generate_anomaly_melt_map(
         keep_year_label_wrapped=True,
         reset_picklefile=False,
         message_below_year="relative to 1990-2020",
-        verbose=True,
     ):
         """Generate a cumulative annual anomaly melt map compared to the baseline climatology period.
 
@@ -1864,7 +1856,7 @@ def generate_anomaly_melt_map(
 
             if mmdd_of_year is None:
                 # Just get the annual anomlay map for that year.
-                anomaly_data = read_annual_melt_anomaly_tif(year=year, verbose=verbose)
+                anomaly_data = read_annual_melt_anomaly_tif(year=year)
             else:
                 datetime_this_year = datetime.datetime(
                     year=year
@@ -1875,7 +1867,6 @@ def generate_anomaly_melt_map(
                 anomaly_data = create_partial_year_melt_anomaly_tif(
                     current_datetime=datetime_this_year,
                     gap_filled=False,
-                    verbose=verbose,
                 )
 
             if anomaly_data is None:
@@ -1966,7 +1957,6 @@ def generate_latest_partial_anomaly_melt_map(
         keep_year_label_wrapped=True,
         reset_picklefile=False,
         message_below_year=None,
-        verbose=True,
     ):
         """Same as generate_anomaly_melt_map, but do it for only a partial year,
         up until the last day of data that we have in the melt array.
@@ -1998,7 +1988,6 @@ def generate_latest_partial_anomaly_melt_map(
             keep_year_label_wrapped=keep_year_label_wrapped,
             reset_picklefile=reset_picklefile,
             message_below_year=message_below_year,
-            verbose=verbose,
         )
 
 
@@ -2014,7 +2003,7 @@ def SPECIAL_make_map_with_borders(year=2020):
         DATA_QGIS_DIR / "basins " / "Antarctic_Regions_v2_interior_borders.shp"
     )
 
-    at = AT_map_generator(fill_pole_hole=False, verbose=True)
+    at = AT_map_generator(fill_pole_hole=False)
     for fmt in ("png", "svg"):
         # for fmt in ("png",):
         fname = os.path.join(
diff --git a/antarctica_today/generate_daily_melt_file.py b/antarctica_today/generate_daily_melt_file.py
index 51d82fd..5e60a96 100644
--- a/antarctica_today/generate_daily_melt_file.py
+++ b/antarctica_today/generate_daily_melt_file.py
@@ -21,6 +21,7 @@
 
 import numpy
 import xarray
+from loguru import logger
 
 from antarctica_today import tb_file_data, write_NSIDC_bin_to_gtif
 from antarctica_today.melt_array_picklefile import get_ice_mask_array
@@ -91,11 +92,13 @@ def generate_new_daily_melt_files(
         # Make sure there's at least one of each file (i.e. exactly one). If not, just skip & continue
         if len(nsidc_fps) == 0:
             if warn_if_missing_files:
-                warnings.warn(
+                msg = (
                     "Warning: At least one NSIDC Tb file on date '"
                     + dt.strftime("%Y%m%d")
                     + "' is missing. Skipping that date."
                 )
+                logger.warning(msg)
+                warnings.warn(msg)
             continue
 
         threshold_file = get_correct_threshold_file(dt)
@@ -119,7 +122,6 @@ def create_daily_melt_file(
     output_bin_filename,
     output_gtif_filename=None,
     Tb_nodata_value=-999,
-    verbose=True,
 ) -> numpy.ndarray:
     """Read input files and generate a daily melt file. Primary function."""
     output_array = read_files_and_generate_melt_array(
@@ -131,7 +133,10 @@ def create_daily_melt_file(
     # Write the output .bin file
     # write_flat_binary.write_array_to_binary(
     write_array_to_binary(
-        output_array, output_bin_filename, numbytes=2, signed=True, verbose=verbose
+        output_array,
+        output_bin_filename,
+        numbytes=2,
+        signed=True,
     )
 
     # Write the output.tif file, if called for
@@ -142,7 +147,6 @@ def create_daily_melt_file(
             resolution=25,
             hemisphere="S",
             nodata=None,
-            verbose=verbose,
         )
 
     return output_array
diff --git a/antarctica_today/generate_plots_for_given_day.py b/antarctica_today/generate_plots_for_given_day.py
index d926498..752f2d9 100644
--- a/antarctica_today/generate_plots_for_given_day.py
+++ b/antarctica_today/generate_plots_for_given_day.py
@@ -106,7 +106,6 @@ def generate_maps_and_plots_for_a_date(
             gap_filled=True,
             dpi=dpi,
             outfile=lineplot_outfile,
-            verbose=True,
         )
 
         # Close the current plots open in matplotlib. (Keeps them from accumulating.)
diff --git a/antarctica_today/progress_bar.py b/antarctica_today/progress_bar.py
index 048e4b3..ffee5f9 100755
--- a/antarctica_today/progress_bar.py
+++ b/antarctica_today/progress_bar.py
@@ -34,10 +34,10 @@ def ProgressBar(
     )
     filledLength = int((length * iteration) // total)
     bar = fill * filledLength + "-" * (length - filledLength)
-    print(f"\r{prefix} |{bar}| {percent}% {suffix}", end=printEnd)
+    print(f"\r{prefix} |{bar}| {percent}% {suffix}", end=printEnd)  # noqa: T201
     # Print New Line on Complete
     if iteration == total:
-        print()
+        print()  # noqa: T201
 
 
 # Sample Usage
diff --git a/antarctica_today/src_baseline/download_NSIDC_Tb_files.py b/antarctica_today/src_baseline/download_NSIDC_Tb_files.py
index 2fecd24..d31d7c0 100644
--- a/antarctica_today/src_baseline/download_NSIDC_Tb_files.py
+++ b/antarctica_today/src_baseline/download_NSIDC_Tb_files.py
@@ -36,6 +36,8 @@
 from getpass import getpass
 from typing import List
 
+from loguru import logger
+
 try:
     from urllib.error import HTTPError, URLError
     from urllib.parse import urlparse
@@ -107,7 +109,7 @@ def get_credentials(url):
         errprefix = "netrc error: "
     except Exception as e:
         if not ("No such file" in str(e)):
-            print("netrc error: {0}".format(str(e)))
+            logger.error("netrc error: {0}".format(str(e)))
         username = None
         password = None
 
@@ -125,7 +127,7 @@ def get_credentials(url):
                 opener = build_opener(HTTPCookieProcessor())
                 opener.open(req)
             except HTTPError:
-                print(errprefix + "Incorrect username or password")
+                logger.error(errprefix + "Incorrect username or password")
                 errprefix = ""
                 credentials = None
                 username = None
@@ -137,7 +139,7 @@ def get_credentials(url):
 def build_version_query_params(version):
     desired_pad_length = 3
     if len(version) > desired_pad_length:
-        print('Version string too long: "{0}"'.format(version))
+        logger.error('Version string too long: "{0}"'.format(version))
         quit()
 
     version = str(int(version))  # Strip off any leading zeros
@@ -178,7 +180,7 @@ def cmr_download(urls, output_dir=None, credentials=None):
         return
 
     url_count = len(urls)
-    print("Downloading {0} files...".format(url_count))
+    logger.info("Downloading {0} files...".format(url_count))
     # credentials = None
 
     for index, url in enumerate(urls, start=1):
@@ -188,7 +190,7 @@ def cmr_download(urls, output_dir=None, credentials=None):
         filename = url.split("/")[-1]
         if output_dir != None:
             filename = os.path.join(output_dir, filename)
-        print(
+        logger.info(
             "{0}/{1}: {2}".format(
                 str(index).zfill(len(str(url_count))), url_count, filename
             )
@@ -205,9 +207,9 @@ def cmr_download(urls, output_dir=None, credentials=None):
             data = opener.open(req).read()
             open(filename, "wb").write(data)
         except HTTPError as e:
-            print("HTTP error {0}, {1}".format(e.code, e.reason))
+            logger.info("HTTP error {0}, {1}".format(e.code, e.reason))
         except URLError as e:
-            print("URL error: {0}".format(e.reason))
+            logger.info("URL error: {0}".format(e.reason))
         except IOError:
             raise
         except KeyboardInterrupt:
@@ -272,7 +274,7 @@ def cmr_search(
         polygon=polygon,
         filename_filter=filename_filter,
     )
-    print("Querying for data:\n\t{0}\n".format(cmr_query_url))
+    logger.info("Querying for data:\n\t{0}\n".format(cmr_query_url))
 
     cmr_scroll_id = None
     ctx = ssl.create_default_context()
@@ -292,21 +294,21 @@ def cmr_search(
                 cmr_scroll_id = headers["cmr-scroll-id"]
                 hits = int(headers["cmr-hits"])
                 if hits > 0:
-                    print("Found {0} matches.".format(hits))
+                    logger.info("Found {0} matches.".format(hits))
                 else:
-                    print("Found no matches.")
+                    logger.info("Found no matches.")
             search_page = response.read()
             search_page = json.loads(search_page.decode("utf-8"))
             url_scroll_results = cmr_filter_urls(search_page)
             if not url_scroll_results:
                 break
             if hits > CMR_PAGE_SIZE:
-                print(".", end="")
+                print(".", end="")  # noqa: T201
                 sys.stdout.flush()
             urls += url_scroll_results
 
         if hits > CMR_PAGE_SIZE:
-            print()
+            print()  # noqa: T201
         return urls
     except KeyboardInterrupt:
         quit()
diff --git a/antarctica_today/src_baseline/extrapolate_thermap_readings_off_map_edge.py b/antarctica_today/src_baseline/extrapolate_thermap_readings_off_map_edge.py
index e4e5b6b..80d1bf4 100644
--- a/antarctica_today/src_baseline/extrapolate_thermap_readings_off_map_edge.py
+++ b/antarctica_today/src_baseline/extrapolate_thermap_readings_off_map_edge.py
@@ -14,6 +14,7 @@
 import os
 
 import numpy
+from loguru import logger
 from osgeo import gdal
 
 ice_mask_tif = "F:/Research/DATA/Antarctica_Today/baseline_datasets/ice_mask.tif"
@@ -38,7 +39,7 @@
 out_ndv = 0.0
 out_array[out_array == tm_ndv] = out_ndv
 
-print(numpy.where(numpy.logical_and((im_array == 1), (tm_array == tm_ndv))))
+logger.info(numpy.where(numpy.logical_and((im_array == 1), (tm_array == tm_ndv))))
 
 # These are the hand-selected pixel values, eight lines total.
 # Five going vertically along Queen Maud Land, extrapolating 1-2 pixels
@@ -85,9 +86,9 @@ def f(x, a, b, c):
     p = numpy.polyfit(known_x, tm_array[known_i, known_j], 2)
     extrapolated_values = f(numpy.array(interp_x), *p)
 
-    print("\n", known_i, interp_i, known_j, interp_j)
-    print(tm_array[known_i, known_j], extrapolated_values)
-    print(known_x, interp_x)
+    logger.info(f"\n {known_i} {interp_i} {known_j} {interp_j}")
+    logger.info(f"{tm_array[known_i, known_j]} {extrapolated_values}")
+    logger.info(f"{known_x} {interp_x}")
 
     # Fill in missing values with extrapolated values
     out_array[interp_i, interp_j] = extrapolated_values
@@ -119,4 +120,4 @@ def f(x, a, b, c):
 ds_out.FlushCache()
 band_out = None
 ds_out = None
-print("\n", thermap_tif_out, "written.")
+logger.info(f"\n {thermap_tif_out} written.")
diff --git a/antarctica_today/src_baseline/plot_thermap_lapse_rate.py b/antarctica_today/src_baseline/plot_thermap_lapse_rate.py
index 5b09ee9..1bde444 100644
--- a/antarctica_today/src_baseline/plot_thermap_lapse_rate.py
+++ b/antarctica_today/src_baseline/plot_thermap_lapse_rate.py
@@ -3,6 +3,7 @@
 import numpy
 import pandas as pd
 import statsmodels.api as sm
+from loguru import logger
 from matplotlib import pyplot as plt
 from matplotlib.axes import Axes
 from matplotlib.figure import Figure
@@ -30,7 +31,7 @@
 
 X = sm.add_constant(elevs)
 model_elev_only = sm.OLS(thermap_df[["Temp"]], X).fit()
-print(model_elev_only.summary())
+logger.info(model_elev_only.summary())
 coefs = model_elev_only.params
 
 
@@ -75,11 +76,11 @@ def plus_minus_op(x):
 X = thermap_df[["REMA_or_Thermap_Elev", "Lat(S)"]]
 Y = thermap_df[["Temp"]]
 
-print("\n=== Statsmodels ===")
+logger.info("=== Statsmodels ===")
 X = sm.add_constant(X)
 model = sm.OLS(Y, X).fit()
 
-print(model.summary())
+logger.info(model.summary())
 coefs = model.params
 
 temps_lat_corrected_75 = temps - coefs["Lat(S)"] * (75 + lats)
@@ -87,7 +88,7 @@ def plus_minus_op(x):
 
 # # Compute a quadratic curve through this line.
 # poly_coefs = numpy.polyfit(elevs, temps_lat_corrected_75, deg=2)
-# print(poly_coefs)
+# logger.info(poly_coefs)
 # # Quadratic trend-line
 # trend_x = numpy.linspace(*min_max_elev, 100)
 # trend_y = poly_coefs[0]*(trend_x**2) + poly_coefs[1]*trend_x + poly_coefs[2]
diff --git a/antarctica_today/src_baseline/resize_25m_grid_by_1000.py b/antarctica_today/src_baseline/resize_25m_grid_by_1000.py
index d9d7e73..fb6b894 100644
--- a/antarctica_today/src_baseline/resize_25m_grid_by_1000.py
+++ b/antarctica_today/src_baseline/resize_25m_grid_by_1000.py
@@ -9,6 +9,7 @@
 
 import os
 
+from loguru import logger
 from osgeo import gdal
 
 infile = "C:/Users/mmacferrin/Dropbox/Research/Antarctica_Today/Dan Dixon/derived/polar_grid_10m_temps_25m_OFF_BY_1000.tif"
@@ -45,4 +46,4 @@
 ds_out.FlushCache()
 band_out = None
 ds_out = None
-print(outfile, "written.")
+logger.info(f"Wrote {outfile}")
diff --git a/antarctica_today/src_baseline/resize_grid_to_reference.py b/antarctica_today/src_baseline/resize_grid_to_reference.py
index 0f7fee1..5ca1063 100644
--- a/antarctica_today/src_baseline/resize_grid_to_reference.py
+++ b/antarctica_today/src_baseline/resize_grid_to_reference.py
@@ -9,24 +9,23 @@
 
 import gdal
 import numpy
+from loguru import logger
 
 
-def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out, verbose=False):
+def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out):
     """I have RACMO & REMA files written out the same grid format & resolution
     as the NSIDC's nsidc-0001 and nsidc-0080 files. But the grid sizes are different
     with different boundaries. This takes a .tif GeoTiff, and a reference Tb GeoTiff,
     and creates a copy of the gtif_in data with the same array size as the gtif_Tb_reference,
     and spits it out to gtif_out. Extra values are filled in with the gtif_in NoDataValue.
     """
-    if verbose:
-        print("Reading", os.path.split(gtif_in)[1])
+    logger.debug(f"Reading {os.path.split(gtif_in)[1]}")
 
     ds_in = gdal.Open(gtif_in, gdal.GA_ReadOnly)
     if ds_in is None:
         raise FileNotFoundError("Gdal could not read input file '{0}'".format(gtif_in))
 
-    if verbose:
-        print("Reading", os.path.split(gtif_reference)[1])
+    logger.debug(f"Reading {os.path.split(gtif_reference)[1]}")
 
     ds_ref = gdal.Open(gtif_reference, gdal.GA_ReadOnly)
     if ds_ref is None:
@@ -66,7 +65,9 @@ def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out, verbose=Fals
         (x_UL_in % x_res_in) == (x_UL_ref % x_res_ref)
         and (y_UL_in % y_res_in) == (y_UL_ref % y_res_ref)
     ):
-        print(
+        msg = "Input grids are not geographically aligned."
+        logger.error(msg)
+        logger.error(
             "X: {0} % {1} = {2}, {3} % {4} = {5}".format(
                 x_UL_in,
                 x_res_in,
@@ -76,7 +77,7 @@ def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out, verbose=Fals
                 x_UL_ref % x_res_ref,
             )
         )
-        print(
+        logger.error(
             "Y: {0} % {1} = {2}, {3} % {4} = {5}".format(
                 y_UL_in,
                 y_res_in,
@@ -86,7 +87,7 @@ def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out, verbose=Fals
                 y_UL_ref % y_res_ref,
             )
         )
-        raise ValueError("Input grids are not geographically aligned.")
+        raise ValueError(msg)
 
     # Create the output array, same shape as the reference array, but same datatype
     # as the source array. Fill with the array_in NDV
@@ -147,8 +148,7 @@ def resize_tif_to_reference_grid(gtif_in, gtif_reference, gtif_out, verbose=Fals
     ds_out.FlushCache()
     ds_out = None
 
-    if verbose:
-        print(os.path.split(gtif_out)[-1], "written.")
+    logger.debug(f"Wrote {os.path.split(gtif_out)[-1]}")
 
     return
 
@@ -160,13 +160,6 @@ def read_and_parse_args():
     parser.add_argument("input_gtif", type=str, help="Source file (.tif)")
     parser.add_argument("reference_gtif", type=str, help="Reference file (.tif)")
     parser.add_argument("output_gtif", type=str, help="Destination file (.tif)")
-    parser.add_argument(
-        "--verbose",
-        "-v",
-        action="store_true",
-        default=False,
-        help="Increase output verbosity.",
-    )
 
     return parser.parse_args()
 
@@ -175,5 +168,7 @@ def read_and_parse_args():
     args = read_and_parse_args()
 
     resize_tif_to_reference_grid(
-        args.input_gtif, args.reference_gtif, args.output_gtif, verbose=args.verbose
+        args.input_gtif,
+        args.reference_gtif,
+        args.output_gtif,
     )
diff --git a/antarctica_today/src_baseline/sample_REMA_elevations.py b/antarctica_today/src_baseline/sample_REMA_elevations.py
index e1d4130..d7bcd41 100644
--- a/antarctica_today/src_baseline/sample_REMA_elevations.py
+++ b/antarctica_today/src_baseline/sample_REMA_elevations.py
@@ -4,6 +4,7 @@
 
 import numpy
 import pandas as pd
+from loguru import logger
 from matplotlib import pyplot as plt
 
 EGM96 = True
@@ -19,7 +20,7 @@
 thermap_df["REMA_Elev"] = [0] * len(thermap_df)
 thermap_df["REMA_or_Thermap_Elev"] = [0] * len(thermap_df)
 
-print(thermap_df.columns)
+logger.info(thermap_df.columns)
 
 elevs = thermap_df["Elev"]
 lons = thermap_df["Lon(W)"]
@@ -40,11 +41,10 @@ def get_rema_elev(row):
         )
     )
     try:
-        print(
+        logger.info(
             "{0:>30s} {1:0.2f} {2:0.2f}, {4:0.2f}*C, {3:0.1f} -> ".format(
                 row["Name"], row["Lat(S)"], row["Lon(W)"], row["Elev"], row["Temp"]
-            ),
-            end="",
+            )
         )
         elev_value = float(return_line.read())
         row["REMA_Elev"] = float(elev_value)
@@ -61,12 +61,12 @@ def get_rema_elev(row):
 
     rema_elevs[idx] = row["REMA_Elev"]  # type: ignore [call-overload]
     rema_or_thermap_elevs[idx] = row["REMA_or_Thermap_Elev"]  # type: ignore [call-overload]
-    print("{0:0.1f}".format(row["REMA_Elev"]))
+    logger.info("-> {0:0.1f}".format(row["REMA_Elev"]))
 
 thermap_df["REMA_Elev"] = rema_elevs
 thermap_df["REMA_or_Thermap_Elev"] = rema_or_thermap_elevs
 
-print("Done")
+logger.info("Done")
 
 fig, ax = plt.subplots(1, 1, figsize=(4, 4))
 ax.set_aspect("equal")
@@ -90,8 +90,8 @@ def get_rema_elev(row):
     )
 
 fig.savefig(fig_outfile, dpi=120)
-print(os.path.split(fig_outfile)[1], "saved.")
+logger.info(f"Saved {os.path.split(fig_outfile)[1]}")
 
 thermap_df.fillna("", inplace=True)
 thermap_df.to_csv(output_csv, index=False, header=True)
-print(os.path.split(output_csv)[1], "saved.")
+logger.info(f"Saved {os.path.split(output_csv)[1]}")
diff --git a/antarctica_today/svgclip.py b/antarctica_today/svgclip.py
index 5d43451..f2605c6 100755
--- a/antarctica_today/svgclip.py
+++ b/antarctica_today/svgclip.py
@@ -32,6 +32,7 @@
 
 import cairo
 import gi
+from loguru import logger
 
 gi.require_version("Rsvg", "2.0")
 gi.require_foreign("cairo")
@@ -58,7 +59,7 @@ def get_bounding_box(svgfile):
 
 def print_info(svgfile):
     bbox = get_bounding_box(svgfile)
-    print(
+    logger.info(
         """
 X:      %f
 Y:      %f
diff --git a/antarctica_today/update_data.py b/antarctica_today/update_data.py
index c478ff8..3bc3126 100644
--- a/antarctica_today/update_data.py
+++ b/antarctica_today/update_data.py
@@ -13,6 +13,7 @@
 import dateutil.parser
 import matplotlib.pyplot
 import numpy
+from loguru import logger
 
 from antarctica_today import (
     generate_daily_melt_file,
@@ -48,7 +49,7 @@ def get_list_of_NSIDC_bin_files_to_import(
         if os.path.splitext(f)[1].lower() == target_extension.strip().lower()
     ]
 
-    # print(len(file_list_all), "total files.")
+    # logger.info(len(file_list_all), "total files.")
 
     # Filter out only the files we want.
     hemisphere_lower = hemisphere.lower()
@@ -63,7 +64,7 @@ def get_list_of_NSIDC_bin_files_to_import(
         + ".bin$"
     )
 
-    # print(search_template)
+    # logger.info(search_template)
     # Create a compiled regular-expression search object
     pattern = re.compile(search_template)
     # Keep only the file names that match the search pattern.
@@ -182,7 +183,7 @@ def update_everything_to_latest_date(
             fnames_this_date[0].suffix.lower() == ".nc"
         ):
             # Read in netCDF file here.
-            # print("Generating melt file {0} from {1}.".format(os.path.basename(melt_bin_fname),
+            # logger.info("Generating melt file {0} from {1}.".format(os.path.basename(melt_bin_fname),
             #                                                   os.path.basename(fnames_this_date[0])
             #                                                   )
             #       )
@@ -261,11 +262,11 @@ def update_everything_to_latest_date(
     new_daily_melt_arrays = []
     new_daily_dts = []
 
-    # print("dt_today:", dt_today)
-    # print("latest_dt_in_array:", latest_dt_in_array)
-    # print(range(1, ((dt_today - latest_dt_in_array).days + 1)))
-    # print(melt_bin_files[-1])
-    # print(melt_bin_paths[-1])
+    # logger.info("dt_today:", dt_today)
+    # logger.info("latest_dt_in_array:", latest_dt_in_array)
+    # logger.info(range(1, ((dt_today - latest_dt_in_array).days + 1)))
+    # logger.info(melt_bin_files[-1])
+    # logger.info(melt_bin_paths[-1])
 
     # For each day, find the .bin file for that day (if it exists) and append it to the list.
     for day_delta in range(1, ((dt_today - latest_dt_in_array).days + 1)):
@@ -288,7 +289,7 @@ def update_everything_to_latest_date(
         daily_melt_array = read_NSIDC_bin_file(
             melt_filepath, element_size=2, return_type=int, signed=True, multiplier=1
         )
-        print(melt_filename, "read.")
+        logger.info(f"Read {melt_filename}")
 
         # Add a 3rd (time) dimension to each array to allow concatenating.
         daily_melt_array.shape = list(daily_melt_array.shape) + [1]
@@ -313,7 +314,7 @@ def update_everything_to_latest_date(
             pickle.dump((melt_array_updated, dt_dict), f)
             f.close()
 
-        print(tb_file_data.model_results_picklefile, "written.")
+        logger.info(f"Wrote {tb_file_data.model_results_picklefile}")
 
     else:
         melt_array_updated = previous_melt_array
@@ -375,7 +376,6 @@ def update_everything_to_latest_date(
             region_num=region_num,
             gap_filled=True,
             outfile=line_plot_outfile,
-            verbose=True,
         )
 
         # Clear the figures just made above in order to not get too many open at once.
@@ -406,7 +406,6 @@ def copy_latest_date_plots_to_date_directory(
     anomaly_maps_dir=map_filedata.anomaly_maps_directory,
     line_plots_dir=tb_file_data.climatology_plots_directory,
     use_symlinks=True,
-    verbose=True,
 ):
     """After running the 'update_everything_to_latest_date()' function, use this to gather all the
     latest-date plots into one location. Put it in a sub-directory of the daily_plots_gathered_dir
@@ -453,8 +452,7 @@ def copy_latest_date_plots_to_date_directory(
     dest_dir_location = os.path.join(dest_parent_dir, date_string)
     if not os.path.exists(dest_dir_location):
         os.mkdir(dest_dir_location)
-        if verbose:
-            print("Created directory '{0}'.".format(dest_dir_location))
+        logger.debug("Created directory '{0}'.".format(dest_dir_location))
 
     for fn in files_to_move:
         src = fn
@@ -467,8 +465,7 @@ def copy_latest_date_plots_to_date_directory(
         else:
             shutil.copyfile(src, dst)
 
-        if verbose:
-            print("{0} -> {1}.".format(src, dst))
+        logger.debug("{0} -> {1}.".format(src, dst))
 
 
 if __name__ == "__main__":
diff --git a/antarctica_today/write_NSIDC_bin_to_gtif.py b/antarctica_today/write_NSIDC_bin_to_gtif.py
index 5126c40..4002e6a 100644
--- a/antarctica_today/write_NSIDC_bin_to_gtif.py
+++ b/antarctica_today/write_NSIDC_bin_to_gtif.py
@@ -10,6 +10,7 @@
 from typing import Type
 
 import numpy
+from loguru import logger
 from osgeo import gdal, osr
 
 from antarctica_today.read_NSIDC_bin_file import read_NSIDC_bin_file
@@ -112,7 +113,6 @@ def output_bin_to_gtif(
     header_size=0,
     resolution=None,
     hemisphere=None,
-    verbose=True,
     nodata=0,
     signed=False,
     multiplier="auto",
@@ -139,9 +139,6 @@ def output_bin_to_gtif(
     hemisphere = "N" or "S"
                  If None, the hemisphere is derived from the nsidc-0001 filename.
 
-    verbose = Verbosity of the output. False will run this silently. True will
-              produce feedback to stdout. (default True)
-
     nodata = Nodata value to put in the geotiff. Defaults to 0.0
 
     return_type = The data type of the geotiff raster band. Defaults to float.
@@ -199,7 +196,6 @@ def output_bin_to_gtif(
         resolution=resolution,
         hemisphere=hemisphere,
         nodata=nodata,
-        verbose=verbose,
     )
 
     return
@@ -219,7 +215,11 @@ def get_nsidc_geotransform(hemisphere, resolution):
 
 
 def output_gtif(
-    array, gtif_file, resolution=25, hemisphere="S", nodata=0, verbose=True
+    array,
+    gtif_file,
+    resolution=25,
+    hemisphere="S",
+    nodata=0,
 ):
     """Take an array, output to a geotiff in the NSIDC resolution specified.
 
@@ -288,8 +288,7 @@ def output_gtif(
     ds.FlushCache()
     ds = None
 
-    if verbose:
-        print(gtif_file, "written.")
+    logger.debug(f"Wrote {gtif_file}")
 
     return
 
@@ -360,13 +359,6 @@ def read_and_parse_args():
         default=False,
         help="Read bin as signed data. Default to unsigned.",
     )
-    parser.add_argument(
-        "--verbose",
-        "-v",
-        action="store_true",
-        default=False,
-        help="Increase output verbosity.",
-    )
 
     return parser.parse_args()
 
@@ -433,5 +425,4 @@ def read_and_parse_args():
         nodata=int(args.nodata),
         return_type=out_type,
         multiplier=multiplier,
-        verbose=args.verbose,
     )
diff --git a/antarctica_today/write_flat_binary.py b/antarctica_today/write_flat_binary.py
index 546713a..b58b711 100644
--- a/antarctica_today/write_flat_binary.py
+++ b/antarctica_today/write_flat_binary.py
@@ -11,6 +11,7 @@
 import os
 
 import numpy
+from loguru import logger
 from osgeo import gdal
 
 
@@ -21,7 +22,6 @@ def write_array_to_binary(
     multiplier=1,
     byteorder="little",
     signed=False,
-    verbose=True,
 ):
     if int(numbytes) not in (1, 2, 4, 8):
         raise ValueError("Numbytes must be one of 1,2,4,8.")
@@ -63,8 +63,7 @@ def write_array_to_binary(
 
     f.close()
 
-    if verbose:
-        print(os.path.split(bin_filename)[-1], "written.")
+    logger.debug(f"Wrote {os.path.split(bin_filename)[-1]}")
 
     return bin_filename
 
@@ -77,14 +76,12 @@ def write_gtif_to_binary(
     multiplier=1,
     byteorder="little",
     signed=False,
-    verbose=True,
 ):
     rasterband = int(rasterband)
     if rasterband < 1:
         raise ValueError("Raster band must be an integer greater than or equal to 1.")
 
-    if verbose:
-        print("Reading", os.path.split(gtif_filename)[-1])
+    logger.debug(f"Reading {os.path.split(gtif_filename)[-1]}")
 
     ds = gdal.Open(gtif_filename, gdal.GA_ReadOnly)
     if ds is None:
@@ -110,7 +107,6 @@ def write_gtif_to_binary(
         multiplier=multiplier,
         byteorder=byteorder,
         signed=signed,
-        verbose=verbose,
     )
 
 
@@ -160,13 +156,6 @@ def read_and_parse_args():
         default=False,
         help="Signed data. Defaults to unsigned. Results will be the same if no negative values are in the array.",
     )
-    parser.add_argument(
-        "--verbose",
-        "-v",
-        action="store_true",
-        default=False,
-        help="Increase output verbosity.",
-    )
 
     return parser.parse_args()
 
@@ -180,5 +169,4 @@ def read_and_parse_args():
         numbytes=args.numbytes,
         multiplier=args.multiplier,
         byteorder=args.byteorder,
-        verbose=args.verbose,
     )