Update README and paths for data_figshare

README.md

@@ -12,7 +12,7 @@ This Repository allows:
 - Displaying and downloading thermodynamic solutions including input files, raw data and grids of density, melt fraction,...
 - Computing basic 1-D isostatic balance for a reference case
 
-While you can run these notebooks on Binder, only a part of the data are accessible from there. Full data are localized on FigShare with the same structure as on the github repository. Therefore, we encourage you to download the full dataset on FigShare and to run jupyter-lab locally in order to be able to display everything.
+While you can run these notebooks on Binder (except topo_cont.ipynb that requires too much RAM), only a part of the data are accessible from there. Full data are localized on FigShare with the same structure as on the github repository. Therefore, we encourage you to download the full dataset on FigShare and to run jupyter-lab locally in order to be able to display everything.
 
 ## Content
 
@@ -56,6 +56,9 @@ $ jupyter labextension install \
     ipygany
 ```
 
+You then need to download the data on Figshare following this [link](https://figshare.com/) (the public link will not be available before publication).
+After that, you have to move the `data_figshare` folder into `SM_ESR_isostasy/`.
+
 Finally run the command below to start the Jupyterlab application. It should
 open a new tab in your browser.
 
@@ -117,7 +120,7 @@ The source code and data assets are under the following licenses:
 
     script: MIT. (LICENSE)
     data:
-        thermodyn: CC-BY-4.0 (LICENSE2)
+        thermodyn: CC-BY-4.0 (LICENSE)
         topography: Data from external sources (described in Datasets section)
 
 For a full description, please visit the README and LICENSE files of each package in the corresponding package folders.

isostasy_calib.ipynb

@@ -10,7 +10,10 @@
     "import numpy as np\n",
     "from scripts import *\n",
     "from scipy.interpolate import griddata\n",
-    "from mpl_toolkits.axes_grid1 import make_axes_locatable"
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "path = './data/thermodyn/'\n",
+    "if os.path.exists('./data_figshare/thermodyn/'):\n",
+    "    path='./data_figshare/thermodyn/'"
    ]
   },
   {
@@ -125,23 +128,23 @@
     "                                                     # in the 2-D thermo-mechanical model\n",
     "    geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.tz'\n",
     "    density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.rhoz'\n",
-    "    path                             = './data/geodyn2d/M1/'\n",
+    "    path_used                        = './data/geodyn2d/M1/'\n",
     "    if printOut:\n",
     "        print(\"Density and temperature profiles from 2-D geodynamic model M1 (P-T dependent densities/perplex)\")\n",
     "else:\n",
     "    deltarho_ridge_correction        = 4.05  # kg/m3 # correction to reproduce the elevation difference \n",
     "                                                     # in the 2-D thermo-mechanical model\n",
     "    geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3_M2.tz'\n",
     "    density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3_M2.rhoz'\n",
-    "    path                             = './data/geodyn2d/M2/'\n",
+    "    path_used                        = './data/geodyn2d/M2/'\n",
     "    if printOut:\n",
     "        print(\"Density and temperature profiles from 2-D geodynamic model M2 (temperature dependent only densities)\")\n",
     "\n",
     "lith1 = geodyn1d.lithosphere('ridge_Oc6_5_NoDiffLayer',printOut=False)\n",
     "lith1.update_materials(book_reflith)\n",
-    "lith1.get_imposed_geotherm(path+geotherm_filename,printOut=printOut)\n",
+    "lith1.get_imposed_geotherm(path_used+geotherm_filename,printOut=printOut)\n",
     "lith1.get_pressure_density_profile(printOut=False,\n",
-    "                                   filename=path+density_profile_filename_ridge,\n",
+    "                                   filename=path_used+density_profile_filename_ridge,\n",
     "                                       drho=deltarho_ridge_correction,\n",
     "                                       zdrho=[6500,125000])\n",
     "\n",
@@ -168,7 +171,7 @@
     "            lith2 = geodyn1d.lithosphere(complith,printOut=False)\n",
     "            lith2.update_materials(book_complith)\n",
     "            lith2.get_steady_state_geotherm(dz=dz,printOut=printOut)\n",
-    "            lith2.get_pressure_density_profile(printOut=False,path='./data/thermodyn/'+perplex_name)\n",
+    "            lith2.get_pressure_density_profile(printOut=False,path=path+perplex_name)\n",
     "            pref  = lith2.get_reference_pressure(ref_depth,printOut=False)\n",
     "            # compute average crustal density for information\n",
     "            rhoc = lith2.P_rho_DepthProfile.rho[lith2.P_rho_DepthProfile.Z<35000]\n",
@@ -213,14 +216,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "id": "82a66fcd-ac1c-42b7-bdbf-c27746ddde16",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c7e0582192414a2ebc9bc778b19edace",
+       "model_id": "c921153b724b43deb171129ab4013e2e",
        "version_major": 2,
        "version_minor": 0
       },
@@ -236,7 +239,7 @@
     "# Calibration of the crustal density and degree of depletion of the continentlal lithospheric mantle\n",
     "# and sensitivity on water depth/sea level\n",
     "#\n",
-    "# Plot of the solution (Figures 10b and S13)\n",
+    "# Plot of the solutions (Figures 10b and S13)\n",
     "#\n",
     "\n",
     "%matplotlib widget\n",
@@ -247,22 +250,26 @@
     "\n",
     "sea_level_ref = -400e0    # meters = reference sea level based on average elevation of continents\n",
     "\n",
+    "# Calibration of the crustal density and degree of depletion of the continentlal lithospheric mantle\n",
+    "#\n",
+    "#refdensity_mantle = \"SP2008_speed3_0.55_bl400_waterload\"          # fixed sea level\n",
+    "refdensity_mantle = \"SP2008_speed3_0.55_waterdepth2750_waterload\" # fixed water depth the ridge\n",
+    "\n",
     "\n",
+    "# Sensitivity on water depth/sea level for temperature dependent only densities\n",
+    "#\n",
+    "# Reference mantle density 3300\n",
     "#refdensity_mantle = \"3300_bl400_waterload\"              # fixed sea level\n",
     "#refdensity_mantle = \"3300_waterdepth2750_waterload\"     # fixed water depth the ridge\n",
-    "\n",
+    "# Reference mantle density 3311\n",
     "#refdensity_mantle = \"3311_bl400_waterload\"              # fixed sea level\n",
     "#refdensity_mantle = \"3311_waterdepth2750_waterload\"     # fixed water depth the ridge\n",
     "\n",
-    "#refdensity_mantle = \"SP2008_speed3_0.55_bl400_waterload\"          # fixed sea level\n",
-    "refdensity_mantle = \"SP2008_speed3_0.55_waterdepth2750_waterload\" # fixed water depth the ridge\n",
-    "\n",
-    "\n",
     "\n",
-    "filename = \"rhoc_lithdepletion_\"+str(refdensity_mantle)+\".txt\"\n",
-    "path     = './data/isostasy_calibration/'\n",
+    "filename  = \"rhoc_lithdepletion_\"+str(refdensity_mantle)+\".txt\"\n",
+    "path_used = './data/isostasy_calibration/'\n",
     "\n",
-    "lines = np.loadtxt(path+filename, comments=\"#\", delimiter=\" \", unpack=False)\n",
+    "lines = np.loadtxt(path_used+filename, comments=\"#\", delimiter=\" \", unpack=False)\n",
     "x  = np.asarray(lines[:,0]) # Average depletion buoyancy CLM\n",
     "y  = np.asarray(lines[:,1]) # Reference crustal density\n",
     "dh = np.asarray(lines[:,2])\n",

isostasy_calib2.ipynb

@@ -104,14 +104,14 @@
     "                                         # in the 2-D thermo-mechanical model\n",
     "geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3_M2.tz'\n",
     "density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3_M2.rhoz'\n",
-    "path                             = './data/geodyn2d/M2/'\n",
+    "path_used                        = './data/geodyn2d/M2/'\n",
     "if printOut:\n",
     "    print(\"Density and temperature profiles from 2-D geodynamic model M2 (temperature dependent only densities)\")\n",
     "lith2 = geodyn1d.lithosphere('ridge_Oc6_5_NoDiffLayer',printOut=False)\n",
     "lith2.update_materials(book_reflith)\n",
-    "lith2.get_imposed_geotherm(path+geotherm_filename,printOut=printOut)\n",
+    "lith2.get_imposed_geotherm(path_used+geotherm_filename,printOut=printOut)\n",
     "lith2.get_pressure_density_profile(printOut=False,\n",
-    "                                   filename=path+density_profile_filename_ridge,\n",
+    "                                   filename=path_used+density_profile_filename_ridge,\n",
     "                                       drho=deltarho_ridge_correction,\n",
     "                                      zdrho=[6500,125000])\n",
     "\n",
@@ -168,7 +168,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "0a05ae308aae4f57918ef28e84afb735",
+       "model_id": "094ca40ee5f4478dae8605419bf8810a",
        "version_major": 2,
        "version_minor": 0
       },
@@ -182,7 +182,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "8ee1db53fdea47f2a367c056ec50dc92",
+       "model_id": "9d87d3c878bc40d196648ddf1cfa98fa",
        "version_major": 2,
        "version_minor": 0
       },
@@ -196,7 +196,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "ff45dd4dda35429983091b79766be621",
+       "model_id": "ac64aabcedc14815812af3e1be2be7c5",
        "version_major": 2,
        "version_minor": 0
       },
@@ -406,7 +406,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "da889c9147a44857a720a99d59a77fe3",
+       "model_id": "9f382472d1d448a9bcccd7a24d510d51",
        "version_major": 2,
        "version_minor": 0
       },

isostasy_ref.ipynb

@@ -9,17 +9,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "from scripts import *"
+    "from scripts import *\n",
+    "path = './data/thermodyn/'\n",
+    "if os.path.exists('./data_figshare/thermodyn/'):\n",
+    "    path='./data_figshare/thermodyn/'"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
@@ -126,7 +129,7 @@
     "lith1 = geodyn1d.lithosphere(reflith,printOut=False)\n",
     "lith1.update_materials(book_reflith)\n",
     "lith1.get_steady_state_geotherm(dz=dz,printOut=True)\n",
-    "lith1.get_pressure_density_profile(printOut=False,path='./data/thermodyn/'+perplex_name)\n",
+    "lith1.get_pressure_density_profile(printOut=False,path=path+perplex_name)\n",
     "pref  = lith1.get_reference_pressure(ref_depth,printOut=True)\n",
     "print('Pref = {0:7.5e} Pa with compensation depth of {1:5.1f} km'.format(pref,ref_depth/1e3))\n",
     "# compute average crustal density for information\n",
@@ -141,28 +144,28 @@
     "                                                 # in the 2-D thermo-mechanical model M1\n",
     "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.tz'\n",
     "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.rhoz'\n",
-    "        path                             = './data/geodyn2d/M1/'\n",
+    "        path_used                        = './data/geodyn2d/M1/'\n",
     "    else:\n",
     "        deltarho_ridge_correction        = 4.05 # kg/m3 # correction to reproduce the elevation difference \n",
     "                                                        # in the 2-D thermo-mechanical model M2\n",
     "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3_M2.tz'\n",
     "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3_M2.rhoz'\n",
-    "        path                             = './data/geodyn2d/M2/'\n",
+    "        path_used                        = './data/geodyn2d/M2/'\n",
     "\n",
     "    water_depth_used = water_depth_at_ridge\n",
     "    lith2            = geodyn1d.lithosphere('ridge_Oc6_5_DiffLayer26',printOut=False)\n",
     "    lith2.update_materials(book_complith)\n",
-    "    lith2.get_imposed_geotherm(path+geotherm_filename,printOut=True)\n",
+    "    lith2.get_imposed_geotherm(path_used+geotherm_filename,printOut=True)\n",
     "    lith2.get_pressure_density_profile(printOut=False,\n",
-    "                                       filename=path+density_profile_filename_ridge,\n",
+    "                                       filename=path_used+density_profile_filename_ridge,\n",
     "                                       drho=deltarho_ridge_correction,\n",
     "                                       zdrho=[6500,125000])\n",
     "else:\n",
     "    water_depth_used = 0\n",
     "    lith2            = geodyn1d.lithosphere(complith,printOut=False)\n",
     "    lith2.update_materials(book_complith)\n",
     "    lith2.get_steady_state_geotherm(printOut=True)\n",
-    "    lith2.get_pressure_density_profile(printOut=False,path='./data/thermodyn/'+perplex_name)\n",
+    "    lith2.get_pressure_density_profile(printOut=False,path=path+perplex_name)\n",
     "    # compute average crustal density for information\n",
     "    rhoc = lith2.P_rho_DepthProfile.rho[lith2.P_rho_DepthProfile.Z<35000]\n",
     "    depc = lith2.P_rho_DepthProfile.Z[lith2.P_rho_DepthProfile.Z<35000]\n",
@@ -217,13 +220,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "8a619d40ff944a08b54cd260928d1fab",
+       "model_id": "c27af325c53d4bf587529e4e4fc9d29e",
        "version_major": 2,
        "version_minor": 0
       },
@@ -237,7 +240,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "16f4421e434c48ac8a4766d51acc5747",
+       "model_id": "3a413b6a37a84912a3f753c52e572c8d",
        "version_major": 2,
        "version_minor": 0
       },
@@ -251,7 +254,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "9587e4fe57eb450e8fe88d142f675591",
+       "model_id": "f1adeb376fdc40ec9c99be767afdbedb",
        "version_major": 2,
        "version_minor": 0
       },

scripts/external.py

@@ -371,6 +371,8 @@ def plot_correlation(x,y,title,filename='correlation.pdf',nbins = 20,xlabel='x',
 #--------------------------- READING DATA ---------------------
 
 def load_spreading_rate(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Half Spreading rate
     nc_spreading_rate = netCDF4.Dataset(path+'rate.3.2.nc')
     #print(nc_spreading_rate.variables.keys())
@@ -382,6 +384,8 @@ def load_spreading_rate(path='./data/topography/'):
     return x,y,spreading_rate
 
 def load_seafloor_ages(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Seafloor age
     nc_age = netCDF4.Dataset(path+'age.3.2.nc')
     #print(nc_age.variables.keys())
@@ -393,6 +397,8 @@ def load_seafloor_ages(path='./data/topography/'):
     return x,y,age
 
 def load_strain_rate(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Strain rate
     nc_strain_rate = netCDF4.Dataset(path+'GSRM_strain_2m.grd')
     #print(nc_strain_rate.variables.keys())
@@ -405,6 +411,8 @@ def load_strain_rate(path='./data/topography/'):
     return x,y,strain_rate
 
 def load_etopo(path='./data/topography/',filtered=True,resolution=2,corrected_from_ice=False):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Etopo 1
     collection = False
     if (filtered==True):
@@ -447,6 +455,8 @@ def load_etopo(path='./data/topography/',filtered=True,resolution=2,corrected_fr
     return x,y,elev
 
 def get_shape_etopo(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Etopo 1
     filename  = 'ETOPO1_Bed_g_gmt4_filtered.grd'
     if (os.path.isfile(path+filename)):
@@ -458,6 +468,8 @@ def get_shape_etopo(path='./data/topography/'):
     return nc_etopo1
 
 def load_hotspots(path='./data/topography/',write_id=False,write_grid=False):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Hot spots
     filename = 'dist_closest_hs.grd'
     if write_grid:
@@ -549,6 +561,8 @@ def load_hotspots(path='./data/topography/',write_id=False,write_grid=False):
     return x,y,dist_closest_hs
 
 def load_crustal_thickness(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # Crustal thickness
     filename = 'crust_thickness_GEMMA_2m.grd'
     if (os.path.isfile(path+filename)):
@@ -569,6 +583,8 @@ def load_crustal_thickness(path='./data/topography/'):
     return x,y,crustal_thickness
 
 def load_lithospheric_thickness(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # lithospheric thickness
     filename = 'lith_ave_no_slabs_-180_180_2m.grd'
     if (os.path.isfile(path+filename)):
@@ -590,6 +606,8 @@ def load_lithospheric_thickness(path='./data/topography/'):
     return x,y,lithospheric_thickness
 
 def load_age_lithos(path='./data/topography/'):
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
     # thermal age continent
     filename = 'mant_age_map_-180_180_qgis_2m.grd'
     if (os.path.isfile(path+filename)):
@@ -616,6 +634,8 @@ def define_MOR_pts(path='./data/topography/',selection_name='MOR_pts_all',distan
 
         selection_name =  'MOR_pts_all' or 'MOR_pts_far_from_hs' or 'MOR_pts_close_to_hs'
     '''
+    if os.path.exists('./data_figshare/topography/'):
+        path='./data_figshare/topography/'
 
     filename = './'+selection_name+'_'+str(distance_between_pts_along_ridges)+'km.dat'