diff --git a/doc/examples/notebooks/hydrology.ipynb b/doc/examples/notebooks/hydrology.ipynb
index 1cdbfb13bcf..e32a39a8506 100644
--- a/doc/examples/notebooks/hydrology.ipynb
+++ b/doc/examples/notebooks/hydrology.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "# Hydrology with GRASS\n",
     "\n",
-    "This is a short introduction to common hydrologic workflows in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.script`, the *Python* API for GRASS, and `grass.jupyter`, an experimental *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
+    "This is a short introduction to common hydrologic workflows in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.tools`, the *Python* API for GRASS, and `grass.jupyter`, a *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
     "\n",
     "This interactive notebook is available online thanks to the [https://mybinder.org](Binder) service. To run the select part (called a *cell*), hit `Shift + Enter`.\n"
    ]
@@ -27,10 +27,8 @@
     "# Import Python standard library and IPython packages we need.\n",
     "import subprocess\n",
     "import sys\n",
-    "import csv\n",
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
-    "from collections import defaultdict\n",
     "\n",
     "# Ask GRASS where its Python packages are.\n",
     "sys.path.append(\n",
@@ -38,14 +36,15 @@
     ")\n",
     "\n",
     "# Import the GRASS packages we need.\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
     "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
+    "tools = Tools(session=session)\n",
     "\n",
     "# Set computational region to elevation raster\n",
-    "gs.run_command(\"g.region\", raster=\"elevation\", flags=\"pg\")"
+    "tools.g_region(raster=\"elevation\")"
    ]
   },
   {
@@ -91,29 +90,26 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"r.fill.dir\",\n",
+    "tools.r_fill_dir(\n",
     "    input=\"elevation\",\n",
     "    output=\"elev_fill1\",\n",
     "    direction=\"dir1\",\n",
     "    areas=\"area1\",\n",
     ")\n",
-    "gs.run_command(\n",
-    "    \"r.fill.dir\",\n",
+    "tools.r_fill_dir(\n",
     "    input=\"elev_fill1\",\n",
     "    output=\"elev_fill2\",\n",
     "    direction=\"dir2\",\n",
     "    areas=\"area2\",\n",
     ")\n",
-    "gs.run_command(\n",
-    "    \"r.fill.dir\",\n",
+    "tools.r_fill_dir(\n",
     "    input=\"elev_fill2\",\n",
     "    output=\"elev_fill3\",\n",
     "    direction=\"dir3\",\n",
     "    areas=\"area3\",\n",
     ")\n",
-    "gs.mapcalc(\"depr_bin = if((elevation-elev_fill3) < 0., 1, null())\")\n",
-    "gs.run_command(\"r.colors\", map=\"depr_bin\", color=\"blues\")"
+    "tools.r_mapcalc(expression=\"depr_bin = if((elevation-elev_fill3) < 0., 1, null())\")\n",
+    "tools.r_colors(map=\"depr_bin\", color=\"blues\")"
    ]
   },
   {
@@ -127,7 +123,6 @@
     "# Display the depressions with InteractiveMap to see how they compare to existing waterbodies\n",
     "depr_map = gj.InteractiveMap()\n",
     "depr_map.add_raster(\"depr_bin\")\n",
-    "depr_map.add_layer_control()\n",
     "depr_map.show()"
    ]
   },
@@ -148,8 +143,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"r.watershed\",\n",
+    "tools.r_watershed(\n",
     "    elevation=\"elevation@PERMANENT\",\n",
     "    drainage=\"drainage\",  # Drainage Direction\n",
     "    accumulation=\"flowacc\",  # Flow Accumulation\n",
@@ -159,7 +153,7 @@
     ")\n",
     "\n",
     "# Convert streams raster to vector\n",
-    "gs.run_command(\"r.to.vect\", input=\"streams\", output=\"streams\", type=\"line\")"
+    "tools.r_to_vect(input=\"streams\", output=\"streams\", type=\"line\")"
    ]
   },
   {
@@ -183,7 +177,7 @@
     "# Note that if the raster is located in a different mapset (for example,\n",
     "# elevation is in PERMANENT, not user1), the `r.colors` will not change\n",
     "# the color in InteractiveMap.\n",
-    "gs.run_command(\"r.colors\", map=\"drainage\", color=\"aspect\")\n",
+    "tools.r_colors(map=\"drainage\", color=\"aspect\")\n",
     "\n",
     "# Add elements to map\n",
     "# We set opacity to 1.0 (default is 0.8) so layers won't interfere with each other.\n",
@@ -193,7 +187,6 @@
     "hydro_map.add_raster(\"watersheds\", opacity=1.0)\n",
     "hydro_map.add_vector(\"streams\")\n",
     "\n",
-    "hydro_map.add_layer_control()\n",
     "\n",
     "hydro_map.show()"
    ]
@@ -214,22 +207,25 @@
    "outputs": [],
    "source": [
     "# Count cells in each watershed\n",
-    "gs.run_command(\n",
-    "    \"r.stats.zonal\",\n",
+    "tools.r_stats_zonal(\n",
     "    base=\"watersheds\",\n",
     "    cover=\"elevation\",\n",
     "    method=\"count\",\n",
     "    output=\"watersheds_count\",\n",
     ")\n",
     "\n",
-    "# Get projection resolution\n",
-    "proj = gs.parse_command(\"g.region\", flags=\"m\")\n",
+    "# Get resolution\n",
+    "region = tools.g_region(flags=\"m\", format=\"json\")\n",
     "\n",
     "# Multiply N-S resollution by E-W resolution to get cell area\n",
-    "cell_area = float(proj[\"nsres\"]) * float(proj[\"ewres\"])\n",
+    "cell_area = region[\"nsres\"] * region[\"ewres\"]\n",
     "\n",
     "# Calculate watersheds areas and convert from m2 to km2\n",
-    "gs.mapcalc(\"'watershed_area' = float('watersheds_count'*{})/1000000\".format(cell_area))"
+    "tools.r_mapcalc(\n",
+    "    expression=\"'watershed_area' = float('watersheds_count'*{})/1000000\".format(\n",
+    "        cell_area\n",
+    "    )\n",
+    ")"
    ]
   },
   {
@@ -246,7 +242,7 @@
    "outputs": [],
    "source": [
     "# Display a map of watershed areas.\n",
-    "gs.run_command(\"r.colors\", map=\"watershed_area\", color=\"plasma\")\n",
+    "tools.r_colors(map=\"watershed_area\", color=\"plasma\")\n",
     "\n",
     "watershed_map = gj.Map()\n",
     "watershed_map.d_rast(map=\"watershed_area\")\n",
@@ -282,7 +278,7 @@
    "outputs": [],
    "source": [
     "# Compute Slope\n",
-    "gs.run_command(\"r.slope.aspect\", elevation=\"elevation\", slope=\"slope\")"
+    "tools.r_slope_aspect(elevation=\"elevation\", slope=\"slope\")"
    ]
   },
   {
@@ -318,22 +314,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "separator = \"|\"\n",
-    "\n",
-    "columns = defaultdict(list)  # each value in each column is appended to a list\n",
-    "\n",
-    "text = gs.read_command(\n",
-    "    \"r.univar\", map=\"elevation\", zones=\"watersheds\", separator=separator, flags=\"t\"\n",
-    ")\n",
-    "reader = csv.DictReader(text.splitlines(), delimiter=separator)\n",
-    "for row in reader:  # read a row as {column1: value1, column2: value2,...}\n",
-    "    for k, v in row.items():  # go over each column name and value\n",
-    "        columns[k].append(v)  # append the value into the appropriate list\n",
-    "        # based on column name k\n",
-    "\n",
-    "watersheds = columns[\"zone\"]\n",
-    "means = np.array(columns[\"mean\"], dtype=np.float32)\n",
-    "stddevs = np.array(columns[\"stddev\"], dtype=np.float32)"
+    "univar = tools.r_univar(map=\"elevation\", zones=\"watersheds\", format=\"json\")\n",
+    "watersheds = [zone[\"zone\"] for zone in univar]\n",
+    "means = [zone[\"mean\"] for zone in univar]\n",
+    "stddevs = [zone[\"stddev\"] for zone in univar]"
    ]
   },
   {
@@ -371,9 +355,7 @@
    "outputs": [],
    "source": [
     "# Convert to vector\n",
-    "gs.run_command(\n",
-    "    \"r.to.vect\", flags=\"s\", input=\"watersheds\", output=\"watersheds_vector\", type=\"area\"\n",
-    ")"
+    "tools.r_to_vect(flags=\"s\", input=\"watersheds\", output=\"watersheds_vector\", type=\"area\")"
    ]
   },
   {
diff --git a/doc/examples/notebooks/jupyter_example.ipynb b/doc/examples/notebooks/jupyter_example.ipynb
index 37bdfa415ac..e35059343a1 100644
--- a/doc/examples/notebooks/jupyter_example.ipynb
+++ b/doc/examples/notebooks/jupyter_example.ipynb
@@ -37,11 +37,12 @@
     ")\n",
     "\n",
     "# Import the GRASS packages we need.\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
-    "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")"
+    "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
+    "tools = Tools(session=session)"
    ]
   },
   {
@@ -60,10 +61,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.parse_command(\"g.region\", raster=\"lakes\", flags=\"pg\")\n",
-    "gs.run_command(\n",
-    "    \"r.buffer\", input=\"lakes\", output=\"lakes_buff\", distances=[60, 120, 240, 500]\n",
-    ")\n",
+    "tools.g_region(raster=\"lakes\", flags=\"pg\")\n",
+    "tools.r_buffer(input=\"lakes\", output=\"lakes_buff\", distances=[60, 120, 240, 500])\n",
     "\n",
     "# Start a Map\n",
     "r_buffer_map = gj.Map()\n",
@@ -93,8 +92,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"v.buffer\", input=\"boundary_state\", output=\"buffer\", distance=-10000)\n",
-    "gs.parse_command(\"g.region\", vector=\"boundary_state\", flags=\"pg\")\n",
+    "tools.v_buffer(input=\"boundary_state\", output=\"buffer\", distance=-10000)\n",
+    "tools.g_region(vector=\"boundary_state\", flags=\"pg\")\n",
     "\n",
     "# Start another Map\n",
     "v_buffer_map = gj.Map()\n",
@@ -134,7 +133,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(gs.read_command(\"g.list\", type=\"all\"))"
+    "print(tools.g_list(type=\"all\").text)"
    ]
   },
   {
@@ -150,7 +149,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(gs.read_command(\"g.search.modules\", flags=\"g\"))"
+    "print(tools.g_search_modules(flags=\"g\").text)"
    ]
   }
  ],
diff --git a/doc/examples/notebooks/jupyter_tutorial.ipynb b/doc/examples/notebooks/jupyter_tutorial.ipynb
index 2d2fe26a7dc..be9d9482b95 100644
--- a/doc/examples/notebooks/jupyter_tutorial.ipynb
+++ b/doc/examples/notebooks/jupyter_tutorial.ipynb
@@ -38,14 +38,17 @@
     ")\n",
     "\n",
     "# Import GRASS packages\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
     "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
     "\n",
+    "# Create a Tools object tied to this session.\n",
+    "tools = Tools(session=session)\n",
+    "\n",
     "# Set computational region to the elevation raster.\n",
-    "gs.run_command(\"g.region\", raster=\"elevation\")"
+    "tools.g_region(raster=\"elevation\")"
    ]
   },
   {
@@ -197,7 +200,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"g.region\", save=\"myregion\", n=224000, s=222000, w=633500, e=637300)\n",
+    "tools.g_region(save=\"myregion\", n=224000, s=222000, w=633500, e=637300)\n",
     "myregion_map = gj.Map(saved_region=\"myregion\")\n",
     "myregion_map.d_rast(map=\"elevation\")\n",
     "myregion_map.d_rast(map=\"lakes\")\n",
@@ -479,7 +482,7 @@
    "source": [
     "series = gj.SeriesMap(height=500)\n",
     "series.add_rasters([\"elevation\", \"elevation_shade\", \"slope\"])\n",
-    "series.add_vectors([\"streams\", \"streets\", \"lakes\"])\n",
+    "series.add_vectors([\"streams\", \"streets\", \"roadsmajor\"])\n",
     "series.d_vect(map=\"streets\")\n",
     "series.d_barscale()\n",
     "series.show()  # Create Slider"
diff --git a/doc/examples/notebooks/parallelization_tutorial.ipynb b/doc/examples/notebooks/parallelization_tutorial.ipynb
index 26c0788dd36..88f2b7638f4 100644
--- a/doc/examples/notebooks/parallelization_tutorial.ipynb
+++ b/doc/examples/notebooks/parallelization_tutorial.ipynb
@@ -30,14 +30,15 @@
     ")\n",
     "\n",
     "# Import GRASS packages\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
     "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
+    "tools = Tools(session=session)\n",
     "\n",
     "# Set computational region to the elevation raster.\n",
-    "gs.run_command(\"g.region\", raster=\"elevation\")"
+    "tools.g_region(raster=\"elevation\")"
    ]
   },
   {
diff --git a/doc/examples/notebooks/scripting_example.ipynb b/doc/examples/notebooks/scripting_example.ipynb
index 178d0258acd..b359cdae393 100644
--- a/doc/examples/notebooks/scripting_example.ipynb
+++ b/doc/examples/notebooks/scripting_example.ipynb
@@ -38,15 +38,17 @@
     "\n",
     "# Import the GRASS packages we need.\n",
     "import grass.script as gs\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Create a GRASS session.\n",
     "session = gs.setup.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
     "\n",
+    "# Create a Tools object tied to this session; overwrite=True lets us re-run cells freely.\n",
+    "tools = Tools(session=session, overwrite=True)\n",
+    "\n",
     "# We want functions to raise exceptions and see standard output of the modules in the notebook.\n",
     "gs.set_raise_on_error(True)\n",
     "gs.set_capture_stderr(True)\n",
-    "# Simply overwrite existing maps like we overwrite Python variable values.\n",
-    "os.environ[\"GRASS_OVERWRITE\"] = \"1\"\n",
     "# Enable map rendering in a notebook.\n",
     "os.environ[\"GRASS_FONT\"] = \"sans\"\n",
     "# Set display modules to render into a file (named map.png by default)\n",
@@ -71,13 +73,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.parse_command(\"g.region\", raster=\"lakes\", flags=\"pg\")\n",
-    "gs.run_command(\n",
-    "    \"r.buffer\", input=\"lakes\", output=\"lakes_buff\", distances=[60, 120, 240, 500]\n",
-    ")\n",
-    "gs.run_command(\"d.erase\")\n",
-    "gs.run_command(\"d.rast\", map=\"lakes_buff\")\n",
-    "gs.run_command(\"d.legend\", raster=\"lakes_buff\", range=(2, 5), at=(80, 100, 2, 10))\n",
+    "tools.g_region(raster=\"lakes\")\n",
+    "tools.r_buffer(input=\"lakes\", output=\"lakes_buff\", distances=[60, 120, 240, 500])\n",
+    "tools.d_erase()\n",
+    "tools.d_rast(map=\"lakes_buff\")\n",
+    "tools.d_legend(raster=\"lakes_buff\", range=(2, 5), at=(80, 100, 2, 10))\n",
     "Image(filename=\"map.png\")"
    ]
   },
@@ -98,17 +98,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"v.buffer\", input=\"boundary_state\", output=\"buffer\", distance=-10000)\n",
-    "gs.parse_command(\"g.region\", vector=\"boundary_state\", flags=\"pg\")\n",
-    "gs.run_command(\"d.erase\")  # erase the display before drawing again\n",
+    "tools.v_buffer(input=\"boundary_state\", output=\"buffer\", distance=-10000)\n",
+    "tools.g_region(vector=\"boundary_state\")\n",
+    "tools.d_erase()  # erase the display before drawing again\n",
     "!rm -f $GRASS_LEGEND_FILE  # and remove the legend file\n",
-    "gs.run_command(\n",
-    "    \"d.vect\", map=\"boundary_state\", fill_color=\"#5A91ED\", legend_label=\"State boundary\"\n",
-    ")\n",
-    "gs.run_command(\n",
-    "    \"d.vect\", map=\"buffer\", fill_color=\"#F8766D\", legend_label=\"Inner portion\"\n",
-    ")\n",
-    "gs.run_command(\"d.legend.vect\", at=(10, 35))\n",
+    "tools.d_vect(map=\"boundary_state\", fill_color=\"#5A91ED\", legend_label=\"State boundary\")\n",
+    "tools.d_vect(map=\"buffer\", fill_color=\"#F8766D\", legend_label=\"Inner portion\")\n",
+    "tools.d_legend_vect(at=(10, 35))\n",
     "Image(filename=\"map.png\")"
    ]
   },
@@ -136,7 +132,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(gs.read_command(\"g.list\", type=\"all\"))"
+    "print(tools.g_list(type=\"all\").text)"
    ]
   },
   {
@@ -152,7 +148,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(gs.read_command(\"g.search.modules\", flags=\"g\"))"
+    "print(tools.g_search_modules(flags=\"g\").text)"
    ]
   },
   {
diff --git a/doc/examples/notebooks/solar_potential.ipynb b/doc/examples/notebooks/solar_potential.ipynb
index f53579955e7..cb2d7b75002 100644
--- a/doc/examples/notebooks/solar_potential.ipynb
+++ b/doc/examples/notebooks/solar_potential.ipynb
@@ -7,7 +7,7 @@
    "source": [
     "# Solar Energy Potential Analysis\n",
     "\n",
-    "This is an introduction to solar radiation modeling in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.script`, the *Python* API for GRASS, and `grass.jupyter`, an experimental *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
+    "This is an introduction to solar radiation modeling in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.tools`, the *Python* API for GRASS, and `grass.jupyter`, a *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
     "\n",
     "This interactive notebook is available online thanks to the [https://mybinder.org](Binder) service. To run the select part (called a *cell*), hit `Shift + Enter`."
    ]
@@ -37,14 +37,15 @@
     ")\n",
     "\n",
     "# Import the GRASS packages we need.\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
     "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
+    "tools = Tools(session=session)\n",
     "\n",
     "# Set computational region to elevation raster\n",
-    "gs.run_command(\"g.region\", raster=\"elevation@PERMANENT\", flags=\"pg\")"
+    "tools.g_region(raster=\"elevation@PERMANENT\")"
    ]
   },
   {
@@ -72,7 +73,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"r.slope.aspect\", elevation=\"elevation\", aspect=\"aspect\", slope=\"slope\")"
+    "tools.r_slope_aspect(elevation=\"elevation\", aspect=\"aspect\", slope=\"slope\")"
    ]
   },
   {
@@ -92,8 +93,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"r.sunmask\",\n",
+    "tools.r_sunmask(\n",
     "    elevation=\"elevation\",\n",
     "    year=\"2001\",\n",
     "    month=\"12\",\n",
@@ -121,8 +121,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"r.sun\",\n",
+    "tools.r_sun(\n",
     "    elevation=\"elevation\",\n",
     "    aspect=\"aspect\",\n",
     "    slope=\"slope\",\n",
@@ -130,8 +129,8 @@
     "    day=\"356\",\n",
     "    time=\"16.25\",\n",
     ")\n",
-    "gs.parse_command(\"r.info\", map=\"incident\", flags=\"g\")\n",
-    "gs.run_command(\"r.colors\", map=\"incident\", co=\"bcyr\", flags=\"e\")"
+    "print(tools.r_info(map=\"incident\").text)\n",
+    "tools.r_colors(map=\"incident\", co=\"bcyr\", flags=\"e\")"
    ]
   },
   {
@@ -163,9 +162,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.mapcalc(\"shadow = if(isnull(incident), 1, null())\")\n",
-    "gs.run_command(\"r.colors\", map=\"shadow\", color=\"grey\")\n",
-    "gs.run_command(\"r.colors\", map=\"elevation\", color=\"elevation\")"
+    "tools.r_mapcalc(expression=\"shadow = if(isnull(incident), 1, null())\")\n",
+    "tools.r_colors(map=\"shadow\", color=\"grey\")\n",
+    "tools.r_colors(map=\"elevation\", color=\"elevation\")"
    ]
   },
   {
@@ -199,8 +198,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"r.sun\",\n",
+    "tools.r_sun(\n",
     "    elevation=\"elevation\",\n",
     "    aspect=\"aspect\",\n",
     "    slope=\"slope\",\n",
@@ -208,10 +206,9 @@
     "    glob_rad=\"winter\",\n",
     "    insol_time=\"its356\",\n",
     ")\n",
-    "gs.run_command(\"r.colors\", map=\"winter\", co=\"gyr\", flags=\"e\")\n",
+    "tools.r_colors(map=\"winter\", co=\"gyr\", flags=\"e\")\n",
     "\n",
-    "gs.run_command(\n",
-    "    \"r.sun\",\n",
+    "tools.r_sun(\n",
     "    elevation=\"elevation\",\n",
     "    aspect=\"aspect\",\n",
     "    slope=\"slope\",\n",
@@ -219,7 +216,7 @@
     "    glob_rad=\"summer\",\n",
     "    insol_time=\"its172\",\n",
     ")\n",
-    "gs.run_command(\"r.colors\", map=\"summer\", co=\"gyr\", flags=\"e\")"
+    "tools.r_colors(map=\"summer\", co=\"gyr\", flags=\"e\")"
    ]
   },
   {
diff --git a/doc/examples/notebooks/temporal.ipynb b/doc/examples/notebooks/temporal.ipynb
index 08423086e92..8016bf28d4f 100644
--- a/doc/examples/notebooks/temporal.ipynb
+++ b/doc/examples/notebooks/temporal.ipynb
@@ -58,11 +58,12 @@
     ")\n",
     "\n",
     "# Import GRASS packages\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
-    "gj.init(\"~/data/nc_basic_spm_grass7/climate_2000_2012\")"
+    "session = gj.init(\"~/data/nc_basic_spm_grass7/climate_2000_2012\")\n",
+    "tools = Tools(session=session)"
    ]
   },
   {
@@ -78,8 +79,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"t.create\",\n",
+    "tools.t_create(\n",
     "    output=\"tempmean\",\n",
     "    type=\"strds\",\n",
     "    temporaltype=\"absolute\",\n",
@@ -87,8 +87,7 @@
     "    description=\"Monthly temperature average in NC [deg C]\",\n",
     ")\n",
     "\n",
-    "gs.run_command(\n",
-    "    \"t.create\",\n",
+    "tools.t_create(\n",
     "    output=\"precip_sum\",\n",
     "    type=\"strds\",\n",
     "    temporaltype=\"absolute\",\n",
@@ -110,13 +109,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "tempmean_list = gs.read_command(\n",
-    "    \"g.list\", type=\"raster\", pattern=\"*tempmean\", separator=\"comma\"\n",
-    ").strip()\n",
+    "tempmean_list = tools.g_list(type=\"raster\", pattern=\"*tempmean\", separator=\"comma\").text\n",
     "\n",
-    "precip_list = gs.read_command(\n",
-    "    \"g.list\", type=\"raster\", pattern=\"*precip\", separator=\"comma\"\n",
-    ").strip()"
+    "precip_list = tools.g_list(type=\"raster\", pattern=\"*precip\", separator=\"comma\").text"
    ]
   },
   {
@@ -132,8 +127,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"t.register\",\n",
+    "tools.t_register(\n",
     "    input=\"tempmean\",\n",
     "    type=\"raster\",\n",
     "    start=\"2000-01-01\",\n",
@@ -142,8 +136,7 @@
     "    flags=\"i\",\n",
     ")\n",
     "\n",
-    "gs.run_command(\n",
-    "    \"t.register\",\n",
+    "tools.t_register(\n",
     "    input=\"precip_sum\",\n",
     "    type=\"raster\",\n",
     "    start=\"2000-01-01\",\n",
@@ -166,8 +159,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"t.rast.extract\",\n",
+    "tools.t_rast_extract(\n",
     "    input=\"precip_sum\",\n",
     "    output=\"precip_sum_2010\",\n",
     "    where=\"start_time >= '2010-01-01' and start_time < '2011-01-01'\",\n",
@@ -187,7 +179,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"t.rast.colors\", input=\"precip_sum_2010\", color=\"precipitation_monthly\")"
+    "tools.t_rast_colors(input=\"precip_sum_2010\", color=\"precipitation_monthly\")"
    ]
   },
   {
@@ -248,15 +240,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"t.unregister\",\n",
+    "tools.t_unregister(\n",
     "    type=\"raster\",\n",
     "    input=\"precip_sum_2010\",\n",
     "    maps=\"2010_02_precip,2010_08_precip\",\n",
     ")\n",
-    "print(\n",
-    "    gs.read_command(\"t.rast.list\", input=\"precip_sum_2010\", columns=\"name,start_time\")\n",
-    ")"
+    "print(tools.t_rast_list(input=\"precip_sum_2010\", columns=\"name,start_time\").text)"
    ]
   },
   {
diff --git a/doc/examples/notebooks/viewshed_analysis.ipynb b/doc/examples/notebooks/viewshed_analysis.ipynb
index 1d7f6e452ff..1702b732805 100644
--- a/doc/examples/notebooks/viewshed_analysis.ipynb
+++ b/doc/examples/notebooks/viewshed_analysis.ipynb
@@ -11,7 +11,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This notebook is a basic introduction to viewshed analysis in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.script`, the *Python* API for GRASS, and `grass.jupyter`, an experimental *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
+    "This notebook is a basic introduction to viewshed analysis in *GRASS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.tools`, the *Python* API for GRASS, and `grass.jupyter`, a *Jupyter Notebook* specific package that helps with the launch of *GRASS* and with displaying maps.\n",
     "\n",
     "This interactive notebook is available online thanks to the [https://mybinder.org](Binder) service. To run the select part (called a *cell*), hit `Shift + Enter`."
    ]
@@ -30,6 +30,7 @@
    "outputs": [],
    "source": [
     "# Import Python standard library and IPython packages we need.\n",
+    "import io\n",
     "import subprocess\n",
     "import sys\n",
     "\n",
@@ -39,14 +40,15 @@
     ")\n",
     "\n",
     "# Import the GRASS packages we need.\n",
-    "import grass.script as gs\n",
     "import grass.jupyter as gj\n",
+    "from grass.tools import Tools\n",
     "\n",
     "# Start GRASS Session\n",
     "session = gj.init(\"~/data/nc_basic_spm_grass7/user1\")\n",
+    "tools = Tools(session=session)\n",
     "\n",
     "# Set computational region to elevation raster\n",
-    "gs.run_command(\"g.region\", raster=\"elevation@PERMANENT\", flags=\"pg\")"
+    "tools.g_region(raster=\"elevation@PERMANENT\")"
    ]
   },
   {
@@ -63,12 +65,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.parse_command(\"g.region\", raster=\"elevation\", flags=\"apg\")\n",
-    "gs.write_command(\n",
-    "    \"v.in.ascii\", input=\"-\", stdin=\"%s|%s\" % (642212, 224767), output=\"viewpoints\"\n",
-    ")\n",
-    "gs.run_command(\n",
-    "    \"r.viewshed\",\n",
+    "tools.v_in_ascii(input=io.StringIO(\"%s|%s\" % (642212, 224767)), output=\"viewpoints\")\n",
+    "tools.r_viewshed(\n",
     "    input=\"elevation\",\n",
     "    output=\"tower_los\",\n",
     "    coordinates=\"642212,224767\",\n",
@@ -94,7 +92,6 @@
     "\n",
     "vs_map.add_raster(\"tower_los\", opacity=0.7)\n",
     "vs_map.add_vector(\"viewpoints\")\n",
-    "vs_map.add_layer_control()\n",
     "\n",
     "vs_map.show()"
    ]
diff --git a/doc/jupyter_intro_images.ipynb b/doc/jupyter_intro_images.ipynb
index 423acbd5e9b..1c3213f0971 100644
--- a/doc/jupyter_intro_images.ipynb
+++ b/doc/jupyter_intro_images.ipynb
@@ -8,7 +8,9 @@
    "outputs": [],
    "source": [
     "import grass.jupyter as gj\n",
-    "import grass.script as gs\n",
+    "from grass.tools import Tools\n",
+    "\n",
+    "tools = Tools(overwrite=True)\n",
     "\n",
     "m = gj.Map(width=500)\n",
     "\n",
@@ -46,7 +48,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"r.relief\", input=\"elevation\", output=\"relief\", overwrite=True)\n",
+    "tools.r_relief(input=\"elevation\", output=\"relief\")\n",
     "m = gj.Map(width=500)\n",
     "\n",
     "# Add a shaded relief map\n",
@@ -167,12 +169,10 @@
    "source": [
     "directions = [0, 90, 180, 270]\n",
     "for azimuth in directions:\n",
-    "    gs.run_command(\n",
-    "        \"r.relief\",\n",
+    "    tools.r_relief(\n",
     "        input=\"elevation\",\n",
     "        output=f\"relief_{azimuth}\",\n",
     "        azimuth=azimuth,\n",
-    "        overwrite=True,\n",
     "        zscale=2,\n",
     "    )\n",
     "maps = [f\"relief_{azimuth}\" for azimuth in directions]"
@@ -215,19 +215,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\"g.region\", n=226040, s=223780, e=639170, w=636190)\n",
-    "gs.run_command(\n",
-    "    \"r.slope.aspect\", elevation=\"elevation\", dx=\"dx\", dy=\"dy\", overwrite=True\n",
-    ")\n",
-    "gs.run_command(\n",
-    "    \"r.sim.water\",\n",
+    "tools.g_region(n=226040, s=223780, e=639170, w=636190)\n",
+    "tools.r_slope_aspect(elevation=\"elevation\", dx=\"dx\", dy=\"dy\")\n",
+    "tools.r_sim_water(\n",
     "    flags=\"t\",\n",
     "    elevation=\"elevation\",\n",
     "    dx=\"dx\",\n",
     "    dy=\"dy\",\n",
     "    depth=\"depth\",\n",
     "    duration=30,\n",
-    "    overwrite=True,\n",
     ")"
    ]
   },
@@ -238,13 +234,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gs.run_command(\n",
-    "    \"t.create\",\n",
+    "tools.t_create(\n",
     "    output=\"depth\",\n",
     "    temporaltype=\"relative\",\n",
     "    title=\"Overland flow depth\",\n",
     "    description=\"Overland flow depth\",\n",
-    "    overwrite=True,\n",
     ")"
    ]
   },
@@ -255,8 +249,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "maps = gs.list_strings(type=\"raster\", pattern=\"depth*\")\n",
-    "gs.run_command(\"t.register\", input=\"depth\", maps=maps, overwrite=True)"
+    "maps = [\n",
+    "    item[\"fullname\"]\n",
+    "    for item in tools.g_list(type=\"raster\", pattern=\"depth*\", format=\"json\")\n",
+    "]\n",
+    "tools.t_register(input=\"depth\", maps=maps)"
    ]
   },
   {