diff --git a/environment.yml b/environment.yml
index 7d8ffed..8f1c8c4 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,6 +8,7 @@ dependencies:
   - requests
   - gdal
   - lastools
+  - geopandas
 # --------- dev dep --------- #
   - pre-commit # hooks for applying linters on commit
   - black # code formatting
diff --git a/pdaltools/add_points_in_las.py b/pdaltools/add_points_in_las.py
new file mode 100644
index 0000000..b47419c
--- /dev/null
+++ b/pdaltools/add_points_in_las.py
@@ -0,0 +1,104 @@
+import argparse
+
+import geopandas
+import numpy as np
+import pdal
+
+from pdaltools.las_info import get_writer_parameters_from_reader_metadata, las_info_metadata, get_bounds_from_header_info
+
+
+def extract_points_from_geo(input_geo: str):
+    file = open(input_geo)
+    df = geopandas.read_file(file)
+    return df.get_coordinates(ignore_index=True, include_z=True)
+
+def point_in_bound(bound_minx, bound_maxx, bound_miny, bound_maxy, pt_x, pt_y):
+    return pt_x >= bound_minx and  pt_x <= bound_maxx and pt_y >= bound_miny and  pt_y <= bound_maxy
+
+def add_points_in_las(input_las: str, input_geo: str, output_las: str, inside_las: bool, values_dimensions: {}):
+    points_geo = extract_points_from_geo(input_geo)
+    pipeline = pdal.Pipeline() | pdal.Reader.las(input_las)
+    pipeline.execute()
+    points_las = pipeline.arrays[0]
+    dimensions = list(points_las.dtype.fields.keys())
+
+    if inside_las:
+        mtd = las_info_metadata(input_las)
+        bound_minx, bound_maxx, bound_miny, bound_maxy = get_bounds_from_header_info(mtd)
+
+    for i in points_geo.index:
+        if inside_las :
+            if not point_in_bound(bound_minx, bound_maxx, bound_miny, bound_maxy, points_geo["x"][i], points_geo["y"][i]):
+                continue
+        pt_las = np.empty(1, dtype=points_las.dtype)
+        pt_las[0][dimensions.index("X")] = points_geo["x"][i]
+        pt_las[0][dimensions.index("Y")] = points_geo["y"][i]
+        pt_las[0][dimensions.index("Z")] = points_geo["z"][i]
+        for val in values_dimensions:
+            pt_las[0][dimensions.index(val)] = values_dimensions[val]
+        points_las = np.append(points_las, pt_las, axis=0)
+
+    params = get_writer_parameters_from_reader_metadata(pipeline.metadata)
+    pipeline_end = pdal.Pipeline(arrays=[points_las])
+    pipeline_end |= pdal.Writer.las(output_las, forward="all", **params)
+    pipeline_end.execute()
+
+
+def parse_args():
+    parser = argparse.ArgumentParser("Add points from geometry file in a las/laz file.")
+    parser.add_argument("--input_file", "-i", type=str, help="Las/Laz input file")
+    parser.add_argument("--output_file", "-o", type=str, help="Las/Laz output file.")
+    parser.add_argument("--input_geo_file", "-g", type=str, help="Geometry input file.")
+    parser.add_argument("--inside_las", "-l", type=str, help="Keep points only inside the las boundary.")
+    parser.add_argument(
+        "--dimensions",
+        "-d",
+        metavar="KEY=VALUE",
+        nargs="+",
+        help="Set a number of key-value pairs corresponding to value "
+        "needed in points added in the output las; key should be included in the input las.",
+    )
+    return parser.parse_args()
+
+
+def is_nature(value, nature):
+    if value is None:
+        return False
+    try:
+        nature(value)
+        return True
+    except:
+        return False
+
+
+def parse_var(s):
+    items = s.split("=")
+    key = items[0].strip()
+    if len(items) > 1:
+        value = "=".join(items[1:])
+        if is_nature(value, int):
+            value = int(value)
+        elif is_nature(value, float):
+            value = float(value)
+    return (key, value)
+
+
+def parse_vars(items):
+    d = {}
+    if items:
+        for item in items:
+            key, value = parse_var(item)
+            d[key] = value
+    return d
+
+
+if __name__ == "__main__":
+    args = parse_args()
+    added_dimensions = parse_vars(args.dimensions)
+    add_points_in_las(
+        input_las=args.input_file,
+        input_geo=args.input_geo_file,
+        output_las=args.input_file if args.output_file is None else args.output_file,
+        inside_las=args.inside_las,
+        values_dimensions=added_dimensions,
+    )
diff --git a/pdaltools/las_remove_dimensions.py b/pdaltools/las_remove_dimensions.py
index ea9e3c3..03931ff 100644
--- a/pdaltools/las_remove_dimensions.py
+++ b/pdaltools/las_remove_dimensions.py
@@ -1,5 +1,4 @@
 import argparse
-import os
 
 import pdal
 from pdaltools.las_info import get_writer_parameters_from_reader_metadata
diff --git a/script/test/test_add_points_in_las.sh b/script/test/test_add_points_in_las.sh
new file mode 100644
index 0000000..02dd90e
--- /dev/null
+++ b/script/test/test_add_points_in_las.sh
@@ -0,0 +1,5 @@
+python -u -m pdaltools.add_points_in_las \
+    --input_file test/data/test_data_77055_627760_LA93_IGN69.LAZ \
+    --output_file test/tmp/addded_cmdline.laz \
+    --input_geo_file test/data/add_points/add_points.geojson \
+    --dimensions "Classification"=64 "Intensity"=1.1
\ No newline at end of file
diff --git a/test/data/add_points/add_points.geojson b/test/data/add_points/add_points.geojson
new file mode 100644
index 0000000..129bf7f
--- /dev/null
+++ b/test/data/add_points/add_points.geojson
@@ -0,0 +1,10 @@
+{
+"type": "FeatureCollection",
+"name": "add_points",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2154" } },
+"features": [
+{ "type": "Feature", "properties": { "id": 1 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770558.858057078556158, 6277569.086137800477445, 0. ], [ 770561.225045961793512, 6277574.660014848224819, 0. ], [ 770566.340796128846705, 6277579.241283654235303, 0. ], [ 770566.340796128846705, 6277585.425996542908251, 0. ], [ 770573.747180699137971, 6277589.167366067878902, 0. ], [ 770582.146173510700464, 6277581.455563577823341, 0. ], [ 770587.032860237522982, 6277586.724022705107927, 0. ], [ 770591.308711123419926, 6277590.541746710427105, 0. ], [ 770591.308711123419926, 6277590.541746710427105, 0. ] ] ] } },
+{ "type": "Feature", "properties": { "id": 2 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770568.402367091737688, 6277561.832462190650403, 0. ], [ 770576.114169582375325, 6277564.123096593655646, 0. ] ] ] } },
+{ "type": "Feature", "properties": { "id": 3 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770587.18556919763796, 6277575.423559648916125, 0. ], [ 770590.926938722841442, 6277581.91369045805186, 0. ], [ 770593.141218645963818, 6277588.938302627764642, 0. ] ] ] } }
+]
+}
diff --git a/test/test_add_points_in_las.py b/test/test_add_points_in_las.py
new file mode 100644
index 0000000..971bef7
--- /dev/null
+++ b/test/test_add_points_in_las.py
@@ -0,0 +1,72 @@
+import pytest
+import os
+import random as rand
+import tempfile
+import math
+
+import pdal
+
+import geopandas as gpd
+from shapely.geometry import Point
+
+from pdaltools import add_points_in_las
+
+numeric_precision = 4
+
+TEST_PATH = os.path.dirname(os.path.abspath(__file__))
+INPUT_DIR = os.path.join(TEST_PATH, "data")
+INPUT_LAS = os.path.join(INPUT_DIR, "test_data_77055_627760_LA93_IGN69.laz")
+
+Xmin = 770575
+Ymin = 6277575
+Zmin = 20
+Size = 20
+
+def distance3D(pt_geo, pt_las):
+    return round(
+        math.sqrt((pt_geo.x - pt_las['X']) ** 2 + (pt_geo.y - pt_las['Y']) ** 2 + (pt_geo.z - pt_las['Z']) ** 2),
+        numeric_precision,
+    )
+
+def add_point_in_las(pt_geo, inside_las):
+    geom = [pt_geo]
+    series = gpd.GeoSeries(geom, crs="2154")
+
+    with tempfile.NamedTemporaryFile(suffix="_geom_tmp.las") as out_las_file:
+        with tempfile.NamedTemporaryFile(suffix="_geom_tmp.geojson") as geom_file:
+            series.to_file(geom_file.name)
+
+            added_dimensions = {"Classification":64, "Intensity":1.}
+            add_points_in_las.add_points_in_las(INPUT_LAS, geom_file.name, out_las_file.name, inside_las, added_dimensions)
+
+            pipeline = pdal.Pipeline() | pdal.Reader.las(out_las_file.name)
+            pipeline.execute()
+            points_las = pipeline.arrays[0]
+            points_las = [e for e in points_las if all(e[val] == added_dimensions[val] for val in added_dimensions)]
+            return points_las
+
+def test_add_point_inside_las():
+    X = Xmin + rand.uniform(0, 1) * Size
+    Y = Ymin + rand.uniform(0, 1) * Size
+    Z = Zmin + rand.uniform(0, 1) * 10
+    pt_geo = Point(X, Y, Z)
+    points_las = add_point_in_las(pt_geo=pt_geo, inside_las=True)
+    assert len(points_las) == 1
+    assert distance3D(pt_geo, points_las[0]) < 1 / numeric_precision
+
+def test_add_point_outside_las_no_control():
+    X = Xmin + rand.uniform(2, 3) * Size
+    Y = Ymin + rand.uniform(0, 1) * Size
+    Z = Zmin + rand.uniform(0, 1) * 10
+    pt_geo = Point(X, Y, Z)
+    points_las = add_point_in_las(pt_geo=pt_geo, inside_las=False)
+    assert len(points_las) == 1
+    assert distance3D(pt_geo, points_las[0]) < 1 / numeric_precision
+
+def test_add_point_outside_las_with_control():
+    X = Xmin + rand.uniform(2, 3) * Size
+    Y = Ymin + rand.uniform(2, 3) * Size
+    Z = Zmin + rand.uniform(0, 1) * 10
+    pt_geo = Point(X, Y, Z)
+    points_las = add_point_in_las(pt_geo=pt_geo, inside_las=True)
+    assert len(points_las) == 0