Merge branch 'jkwang/fireinice_python_plotting' into 'main'

Adds the Python code for creating two of the FireInIce visualizations.

See merge request !70

.gitignore

@@ -68,3 +68,22 @@ beaufortSea/images/*
 # sbtools 
 beaufortSea/.Renviron
 
+
+#fire and ice
+fireInIce/*/out
+fireInIce/*/tmp
+fireInIce/*/src/__pycache__
+
+#remove cache and dataout
+fireInIce/cache
+fireInIce/data_out
+
+#remove snakemake files
+fireInIce/.snakemake
+
+#exceptions Spatial
+!fireInIce/2_visualize/in/*.cpg
+!fireInIce/2_visualize/in/*.dbf
+!fireInIce/2_visualize/in/*.prj
+!fireInIce/2_visualize/in/*.shp
+!fireInIce/2_visualize/in/*.shx
\ No newline at end of file

fireInIce/0_fetch/src/download_dem.py

+import os
+import geopandas as gpd
+from bmi_topography import Topography
+import py3dep
+from shapely.geometry import Polygon
+
+def snakemake_type_exists(snakemake_type,string,default_input):
+    if hasattr(snakemake_type, string):
+        return snakemake_type[string]
+    else:
+        return default_input
+
+def get_shapefile_extent(shapefile):
+    west, south, east, north = shapefile.geometry.total_bounds
+    length = east - west
+    height = north - south
+    return west, east, length, south, north, height
+
+
+def buffer_shapefile(west, east, length, south, north, height, buffer):
+    west -= buffer / 100.0 * length
+    south -= buffer / 100.0 * height
+    east += buffer / 100.0 * length
+    north += buffer / 100.0 * height
+    return west, east, south, north
+
+
+def bmi_download_dem(extent_shpfile, demfile, data_product, dem_product, buffer):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(demfile)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # load shapefile
+    shp = gpd.read_file(extent_shpfile)
+    # get extent of shapefile #assums WGS84 or some sort of geographic coordinate system, Potential FIX: could convert if necessary.
+    west, east, length, south, north, height = get_shapefile_extent(shp)
+    buff_west, buff_east, buff_south, buff_north = buffer_shapefile(
+        west, east, length, south, north, height, buffer
+    )
+
+    # set parameters of dem download
+    params = Topography.DEFAULT.copy()
+    params["data_type"] = data_product
+    params["dem_type"] = dem_product
+    params["output_format"] = "GTiff"
+    params["west"] = buff_west
+    params["south"] = buff_south
+    params["east"] = buff_east
+    params["north"] = buff_north
+    params["cache_dir"] = out_dir
+
+    # set api key if exists (you need this when you run out of free requests)
+    if os.path.isfile("opentopography_api_key.txt"):
+        with open("opentopography_api_key.txt", "r") as f:
+            api_key = f.read()
+        params["api_key"] = api_key
+
+    # setup request and fetch
+    topo_request = Topography(**params)
+    topo_request.fetch()
+
+    # rename the file
+    os.rename(
+        out_dir
+        + "/"
+        + str(params["dem_type"])
+        + "_"
+        + str(params["south"])
+        + "_"
+        + str(params["west"])
+        + "_"
+        + str(params["north"])
+        + "_"
+        + str(params["east"])
+        + ".tif",
+        out_dir + "/dem.tif",
+    )
+
+def py3dep_download_dem(extent_shpfile, demfile, py3dep_resolution, buffer):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(demfile)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # load shapefile
+    shp = gpd.read_file(extent_shpfile)
+    # get extent of shapefile #assums WGS84 or some sort of geographic coordinate system, Potential FIX: could convert if necessary.
+    west, east, length, south, north, height = get_shapefile_extent(shp)
+    buff_west, buff_east, buff_south, buff_north = buffer_shapefile(
+        west, east, length, south, north, height, buffer
+    )
+
+    #specify geometry to download dem from
+    buffered_domain_polygon = Polygon(
+        zip(
+            [buff_west, buff_west, buff_east, buff_east],
+            [buff_north, buff_south, buff_south, buff_north],
+        )
+    )
+
+    #download dem
+    dem = py3dep.get_map("DEM", buffered_domain_polygon, resolution = py3dep_resolution, geo_crs = shp.crs, crs = 4326)
+
+    #save dem
+    dem.rio.to_raster(out_dir + "/dem.tif")
+
+if __name__ == "__main__":
+    data_product = snakemake_type_exists(snakemake.params,"data_product","/API/globaldem")
+    dem_product = snakemake_type_exists(snakemake.params,"dem_product","NASADEM")
+    py3dep_resolution = snakemake_type_exists(snakemake.params,"py3dep_resolution",30)
+    buffer = snakemake.params["buffer"]
+    extent_shpfile = snakemake.input["extent_shpfile"]
+    demfile = snakemake.output["demfile"]
+    download_mode = snakemake_type_exists(snakemake.params,"download_mode","bmi")
+    if download_mode == "bmi":
+        bmi_download_dem(extent_shpfile, demfile, data_product, dem_product, buffer)
+    elif download_mode == "py3dep":
+        py3dep_download_dem(extent_shpfile, demfile, py3dep_resolution, buffer)
+    else:
+        raise Exception("Incorrect download mode. Please use bmi or py3dep.")

fireInIce/0_fetch/src/download_dem_bounding_box.py

+import os
+import numpy as np
+import geopandas as gpd
+from bmi_topography import Topography
+from shapely.geometry import Polygon
+
+def snakemake_type_exists(snakemake_type,string,default_input):
+    if hasattr(snakemake_type, string):
+        return snakemake_type[string]
+    else:
+        return default_input
+
+def bmi_download_dem(UL_corner, LR_corner,extent_shpfile, demfile, data_product, dem_product, buffer):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(demfile)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # set parameters of dem download
+    params = Topography.DEFAULT.copy()
+    params["data_type"] = data_product
+    params["dem_type"] = dem_product
+    params["output_format"] = "GTiff"
+    params["west"] = float(UL_corner[1])
+    params["south"] = float(LR_corner[0])
+    params["east"] = float(LR_corner[1])
+    params["north"] = float(UL_corner[0])
+    params["cache_dir"] = out_dir
+
+    # set api key if exists (you need this when you run out of free requests)
+    if os.path.isfile("opentopography_api_key.txt"):
+        with open("opentopography_api_key.txt", "r") as f:
+            api_key = f.read()
+        params["api_key"] = api_key
+
+    # setup request and fetch
+    topo_request = Topography(**params)
+    topo_request.fetch()
+
+    # rename the file 
+    os.rename(
+        out_dir
+        + "/"
+        + str(params["dem_type"])
+        + "_"
+        + str(params["south"])
+        + "_"
+        + str(params["west"])
+        + "_"
+        + str(params["north"])
+        + "_"
+        + str(params["east"])
+        + ".tif",
+        demfile,
+    )
+
+    domain_geom = Polygon(
+        zip(
+            [UL_corner[1],UL_corner[1], LR_corner[1],LR_corner[1]],
+            [UL_corner[0], LR_corner[0], LR_corner[0], UL_corner[0]],
+        )
+    )
+    domain_polygon = gpd.GeoDataFrame(index=[0], crs="EPSG:4326", geometry=[domain_geom])
+    domain_polygon.to_file(extent_shpfile)
+
+
+if __name__ == "__main__":
+    data_product = snakemake_type_exists(snakemake.params,"data_product","/API/globaldem")
+    dem_product = snakemake_type_exists(snakemake.params,"dem_product","NASADEM")
+    buffer = snakemake.params["buffer"]
+    UL_corner = snakemake.params["UL_corner"]
+    LR_corner = snakemake.params["LR_corner"]
+    extent_shpfile = snakemake.output["extent_shpfile"]
+    demfile = snakemake.output["demfile"]
+    download_mode = snakemake_type_exists(snakemake.params,"download_mode","bmi")
+    bmi_download_dem(UL_corner, LR_corner, extent_shpfile, demfile, data_product, dem_product, buffer)

fireInIce/0_fetch/src/download_sb.py

+import os
+import sciencebasepy
+
+def main(sb_id,zipfile):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(zipfile)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+    # download the file
+    sb = sciencebasepy.SbSession()
+    sb_item = sb.get_item(sb_id)
+    sb.get_item_files_zip(sb_item,out_dir)
+
+if __name__ == "__main__":
+    sb_id = snakemake.params["sb_id"]
+    zipfile = snakemake.output["zipfile"]
+    main(sb_id,zipfile)
\ No newline at end of file

fireInIce/0_fetch/src/get_wbd.py

+import os
+import pandas as pd
+import geopandas as gpd
+from pygeohydro import WBD
+
+def main(huc_list, extent_shpfile):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(extent_shpfile)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    #download huc
+    wbd_list = []
+    for huc in huc_list:
+        huc_level = str(len(huc))
+        wbd = WBD("huc" + huc_level)
+        wbd_list += [wbd.byids("huc" + huc_level, [huc]).dissolve()]
+    
+    extent = gpd.GeoDataFrame(pd.concat(wbd_list)).dissolve()
+    extent.to_file(extent_shpfile)
+
+if __name__ == "__main__":
+    huc_list = snakemake.params["huc_list"]
+    extent_shpfile = snakemake.output["extent_shpfile"]
+    main(huc_list, extent_shpfile)

fireInIce/0_fetch/src/unzip_file.py

+import os
+import zipfile as zf
+import py7zr
+
+def unzip_file(zipfile, out_dir):
+    if zipfile.endswith(".7z"):
+        with py7zr.SevenZipFile(zipfile) as zip_ref:
+            zip_ref.extractall(out_dir)
+    elif zipfile.endswith (".zip"):
+        with zf.ZipFile(zipfile, "r") as zip_ref:
+            zip_ref.extractall(out_dir)
+
+def main(zipfile):
+    # make out_dir if it doesn't exist
+    tmp_dir = os.path.dirname(zipfile)
+    if zipfile.endswith(".7z"):
+        out_dir = tmp_dir + '/' + os.path.basename(zipfile[:-3])
+    elif zipfile.endswith (".zip"):
+        out_dir = tmp_dir + '/' + os.path.basename(zipfile[:-4])
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+    # unzip the file
+    unzip_file(zipfile, out_dir)
+
+
+if __name__ == "__main__":
+    zipfile = snakemake.input["zipfile"]
+    main(zipfile)

fireInIce/1_process/src/calc_bedrock.py

+
+import os
+import numpy as np
+from osgeo import gdal, osr
+
+def calc_bedrock(dem_file,ice_thickness_file,ice_height_file,bedrock_file):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(bedrock_file)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # read in the dem
+    ds = gdal.Open(dem_file)
+    ds_it = gdal.Open(ice_thickness_file)
+
+    #get metadata
+    gt = ds.GetGeoTransform()
+    cs = ds.GetProjection()
+    dtype = ds.GetRasterBand(1).DataType
+    cellsx = ds.GetRasterBand(1).XSize
+    cellsy = ds.GetRasterBand(1).YSize
+    NaN = ds.GetRasterBand(1).GetNoDataValue()
+    if NaN == None:
+        NaN = -9999
+
+    #convert to numpy array
+    dem = np.array(ds.GetRasterBand(1).ReadAsArray())
+    it = np.array(ds_it.GetRasterBand(1).ReadAsArray())
+
+    #calc bedrock
+    bedrock = dem - it
+
+    #ice_height
+    ih = np.zeros_like(dem) - 100000.
+    ih[it!=0] = dem[it!=0]
+    
+    driver = gdal.GetDriverByName('GTiff')
+    dataset = driver.Create(ice_height_file,cellsx,cellsy, 1, gdal.GDT_Float32)
+    dataset.GetRasterBand(1).WriteArray(ih)
+    dataset.GetRasterBand(1).SetNoDataValue(NaN)
+
+    dataset.SetGeoTransform(gt)
+    dataset.SetProjection(cs)
+    dataset.FlushCache()
+    dataset=None
+
+    dataset = driver.Create(bedrock_file,cellsx,cellsy, 1, gdal.GDT_Float32)
+    dataset.GetRasterBand(1).WriteArray(bedrock)
+    dataset.GetRasterBand(1).SetNoDataValue(NaN)
+
+    dataset.SetGeoTransform(gt)
+    dataset.SetProjection(cs)
+    dataset.FlushCache()
+    dataset=None
+
+if __name__ == "__main__":
+    dem_file = snakemake.input["dem_file"]
+    ice_thickness_file = snakemake.input["ice_thickness_file"]
+    ice_height_file = snakemake.output["ice_height_file"]
+    bedrock_file = snakemake.output["bedrock_file"]
+    calc_bedrock(dem_file,ice_thickness_file,ice_height_file,bedrock_file)
\ No newline at end of file

fireInIce/1_process/src/cores_to_shp.py

+import pandas as pd
+import geopandas as gpd
+def main(core_data,core_shpfile):
+    cd = pd.read_excel(core_data)
+    cd.loc[2, 'Latitude'] = "58° 51’ 16.1” N" #datafix 2016, C3
+    cd.loc[2, 'Longitude'] = "134° 10’ 31.1” W" #datafix 2016, C3
+    
+    core_name = []
+    Lat = []
+    Long = []
+    for i in range(0,len(cd)):
+        core_name += [cd['U.S. Geological Survey site name'][i]]
+        Lat += [convert_coor_to_dec(cd['Latitude'][i])]
+        Long += [convert_coor_to_dec(cd['Longitude'][i])]
+        print (cd['U.S. Geological Survey site name'][i],cd['Latitude'][i],cd['Longitude'][i],\
+            convert_coor_to_dec(cd['Latitude'][i]),convert_coor_to_dec(cd['Longitude'][i]))
+
+    df = pd.DataFrame(
+        {
+            "Core Name": core_name,
+            "Latitude": Lat,
+            "Longitude": Long,
+        }
+    )
+
+    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude), crs="EPSG:4326")
+
+    gdf.to_file(core_shpfile, driver='ESRI Shapefile')
+
+
+import re
+def convert_coor_to_dec(coor):
+    deg, minutes, seconds, direction =  re.split('[°’”]', coor)
+    return(float(deg) + float(minutes)/60 + float(seconds)/(60*60)) * (-1 if direction in [' W', ' S', ' W', ' S'] else 1)
+
+if __name__ == "__main__":
+    core_data = snakemake.input["core_data"]
+    core_shpfile = snakemake.output["core_shpfile"]
+    main(
+        core_data,
+        core_shpfile
+    )
\ No newline at end of file

fireInIce/1_process/src/regrid.py

+
+#code modified from https://stackoverflow.com/questions/10454316/how-to-project-and-resample-a-grid-to-match-another-grid-with-gdal-python
+from osgeo import gdal, gdalconst
+import os
+
+def main(src_file,match_file,dst_file):
+    out_dir = os.path.dirname(dst_file)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # Source
+    src = gdal.Open(src_file, gdalconst.GA_ReadOnly)
+    src_proj = src.GetProjection()
+
+    # We want a section of source that matches this:
+    match_ds = gdal.Open(match_file, gdalconst.GA_ReadOnly)
+    match_proj = match_ds.GetProjection()
+    match_geotrans = match_ds.GetGeoTransform()
+    wide = match_ds.RasterXSize
+    high = match_ds.RasterYSize
+
+    # Output / destination
+    dst = gdal.GetDriverByName('GTiff').Create(dst_file, wide, high, 1, gdalconst.GDT_Float32)
+    dst.SetGeoTransform( match_geotrans )
+    dst.SetProjection( match_proj)
+
+    # Do the work
+    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
+
+    del dst # Flush
+
+if __name__ == "__main__":
+    src_file = snakemake.input["src_file"]
+    match_file = snakemake.input["match_file"]
+    dst_file = snakemake.output["dst_file"]
+    main(src_file,match_file,dst_file)
\ No newline at end of file

fireInIce/1_process/src/remove_nans.py

+from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
+import os
+import numpy as np
+from osgeo import gdal, osr
+
+x_nei = [-1,0,1,-1,1,-1,0,1]
+y_nei = [1,1,1,0,0,-1,-1,-1]
+
+def seed_ice_hole(it,cellsx,cellsy,it_threshold):
+    for x in range(0,cellsx):
+        for y in range(0,cellsy):
+            if it[y,x] > it_threshold:
+                for nei in range(0,8):
+                    if it[y+y_nei[nei],x+x_nei[nei]] == 0:
+                        it[y+y_nei[nei],x+x_nei[nei]] = -9999
+    return it
+            
+def fill_ice_hole(it,cellsx,cellsy):
+    trigger = 0
+    fill_iteration = 0
+    while trigger == 0:
+        trigger = 1
+        for x in range(0,cellsx):
+            for y in range(0,cellsy):
+                if it[y,x] == -9999:
+                    trigger = 0
+                    mini_trigger = 0
+                    for nei in range(0,8):
+                        if it[y+y_nei[nei],x+x_nei[nei]] == 0:
+                            mini_trigger = 1
+                            it[y+y_nei[nei],x+x_nei[nei]] = -9999
+                    if mini_trigger == 0:
+                        it[y,x] = -999999
+        fill_iteration += 1
+        print ("filling ice, try " + str(fill_iteration) + "...")
+    return it
+
+def remove_nans(ice_thickness_file,fixed_ice_thickness_file,convolution_window):
+    # make out_dir if it doesn't exist
+    out_dir = os.path.dirname(fixed_ice_thickness_file)
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+
+    # read in the dem
+    ds_it = gdal.Open(ice_thickness_file)
+
+    #get metadata
+    gt = ds_it.GetGeoTransform()
+    cs = ds_it.GetProjection()
+    dtype = ds_it.GetRasterBand(1).DataType
+    cellsx = ds_it.GetRasterBand(1).XSize
+    cellsy = ds_it.GetRasterBand(1).YSize
+    NaN =0
+
+    #convert to numpy array
+    it = np.array(ds_it.GetRasterBand(1).ReadAsArray())
+    
+    #convert zero to NaN
+    it = seed_ice_hole(it,cellsx,cellsy,750.)
+    it = fill_ice_hole(it,cellsx,cellsy)
+    it[it==-999999] = np.nan
+
+    #set kernal to one neighbor
+    kernel = Gaussian2DKernel(x_stddev=convolution_window)
+
+    # create a "fixed" image with NaNs replaced by interpolated values
+    fixed_it= interpolate_replace_nans(it, kernel)
+
+    driver = gdal.GetDriverByName('GTiff')
+    dataset = driver.Create(fixed_ice_thickness_file,cellsx,cellsy, 1, gdal.GDT_Float32)
+    dataset.GetRasterBand(1).WriteArray(fixed_it)
+    dataset.GetRasterBand(1).SetNoDataValue(NaN)
+    dataset.SetGeoTransform(gt)
+    dataset.SetProjection(cs)
+    dataset.FlushCache()
+    dataset=None
+
+if __name__ == "__main__":
+    convolution_window = snakemake.params["convolution_window"]
+    ice_thickness_file = snakemake.input["ice_thickness_file"]
+    fixed_ice_thickness_file = snakemake.output["ice_thickness_file"]
+    remove_nans(ice_thickness_file,fixed_ice_thickness_file,convolution_window)
\ No newline at end of file

fireInIce/2_visualize/in/extent_regional_fires.cpg

+UTF-8
\ No newline at end of file

fireInIce/2_visualize/in/extent_regional_fires.prj

+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
\ No newline at end of file

fireInIce/2_visualize/in/photo_locations.cpg

+UTF-8
\ No newline at end of file

fireInIce/2_visualize/in/photo_locations.prj

+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
\ No newline at end of file