diff --git a/README.md b/README.md
index 7f35995..0c17637 100644
--- a/README.md
+++ b/README.md
@@ -112,6 +112,41 @@ Minimalistic C++ app example: [https://github.com/adonmo/meos-cpp-example](https
 
 C++ API Reference: [https://adonmo.github.io/meos/](https://adonmo.github.io/meos/)
 
+## C
+
+MEOS can also be accessed from C.
+
+Here is an example code on how to do this:
+
+```c
+#include "meos/meos_c.h"
+#include <stdio.h>
+
+int main(int argc, char* argv[]) {
+    struct MEOS_TBox* tbox = MEOS_newTBox(12.34, 56.78);
+    double xmin = MEOS_TBox_xmin(tbox);
+    printf("%lf\n", xmin);
+    MEOS_deleteTBox(tbox);
+}
+```
+To compile the above code, you would need `libmeos_c.so`. If you do not have it, build it first by running:
+```sh
+cmake -B build/capi -S capi
+cmake --build build/capi
+```
+
+After saving the above file as `main.c`, run the following commands to build and run it:
+
+```sh
+gcc -I/<path_to_meos>/capi/include -c main.c -o main.o
+g++ -L/<path_to_meos>/build/capi main.o -l:libmeos_c.so -o main
+LD_LIBRARY_PATH=/<path_to_meos>/build/capi ./main
+```
+If you see the following output, you have succesfully used MEOS from C!
+```
+12.340000
+```
+
 ## Contributing
 
 Issues and pull requests are welcome.
diff --git a/capi/CMakeLists.txt b/capi/CMakeLists.txt
new file mode 100644
index 0000000..9e77d8e
--- /dev/null
+++ b/capi/CMakeLists.txt
@@ -0,0 +1,52 @@
+cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
+
+project(
+  meos_c
+  LANGUAGES C CXX
+)
+# Get version from VERSION file instead of hardcoding here
+# This is so that we can maintain a single version across C/C++ and Python
+file(STRINGS ../VERSION LIBMEOS_VERSION)
+set(PROJECT_VERSION ${LIBMEOS_VERSION})
+message(STATUS "MEOS: v${LIBMEOS_VERSION}")
+
+# ---- Dependencies ----
+
+include(../cmake/CPM.cmake)
+
+CPMAddPackage(
+  NAME libmeos
+  SOURCE_DIR ${CMAKE_CURRENT_LIST_DIR}/..
+)
+
+# ---- Create library ----
+
+file(GLOB_RECURSE headers CONFIGURE_DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/include/*.h")
+file(GLOB_RECURSE sources CONFIGURE_DEPENDS "${CMAKE_CURRENT_SOURCE_DIR}/source/*.cpp")
+add_library(meos_c SHARED ${headers} ${sources})
+target_link_libraries(meos_c libmeos)
+
+if(MSVC)
+    add_definitions(-D_CRT_SECURE_NO_WARNINGS)
+endif()
+
+# being a cross-platform target, we enforce standards conformance on MSVC
+target_compile_options(meos_c PUBLIC "$<$<BOOL:${MSVC}>:/permissive->")
+
+target_include_directories(meos_c
+  PUBLIC
+    $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:include/${PROJECT_NAME}-${PROJECT_VERSION}>
+)
+
+message(${PROJECT_SOURCE_DIR}/include)
+
+packageProject(
+  NAME ${PROJECT_NAME}
+  VERSION ${PROJECT_VERSION}
+  BINARY_DIR ${PROJECT_BINARY_DIR}
+  INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include
+  INCLUDE_DESTINATION include/${PROJECT_NAME}-${PROJECT_VERSION}
+  VERSION_HEADER "${VERSION_HEADER_LOCATION}"
+  DEPENDENCIES ""
+)
diff --git a/capi/include/meos/meos_c.h b/capi/include/meos/meos_c.h
new file mode 100644
index 0000000..46f65ff
--- /dev/null
+++ b/capi/include/meos/meos_c.h
@@ -0,0 +1,104 @@
+#pragma once
+
+#include <string.h>
+#include <time.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct MEOS_Period MEOS_Period;
+typedef struct MEOS_TBox MEOS_TBox;
+typedef struct MEOS_TFloatInst MEOS_TFloatInst;
+typedef struct MEOS_TFloatSeq MEOS_TFloatSeq;
+typedef struct MEOS_TFloatSeqSet MEOS_TFloatSeqSet;
+
+typedef enum {
+  MEOS_Interpolation_Stepwise,
+  MEOS_Interpolation_Linear,
+} MEOS_Interpolation;
+
+/*****************************************************************************
+ * TBox
+ *****************************************************************************/
+
+MEOS_TBox *MEOS_newTBox(double const xmin, double const xmax);
+
+double MEOS_TBox_xmin(MEOS_TBox *tbox);
+double MEOS_TBox_xmax(MEOS_TBox *tbox);
+
+void MEOS_deleteTBox(MEOS_TBox *tbox);
+
+/*****************************************************************************
+ * Period
+ *****************************************************************************/
+
+// Constructor with a string
+MEOS_Period *MEOS_newPeriod(char *serialized);
+// Constructor with two timestamps
+MEOS_Period *MEOS_newPeriod_TT(time_t lower, time_t upper);
+// Constructor with two timestamps and two boolean bounds
+MEOS_Period *MEOS_newPeriod_TTBB(time_t lower, time_t upper, bool lower_inc, bool upper_inc);
+
+char *MEOS_Period_str(MEOS_Period *period);
+
+void MEOS_deletePeriod(MEOS_Period *period);
+
+/*****************************************************************************
+ * TFloatInst
+ *****************************************************************************/
+
+// Constructor with a string
+MEOS_TFloatInst *MEOS_newTFloatInst(char *serialized);
+// Constructor with a value and timestamp
+MEOS_TFloatInst *MEOS_newTFloatInst_VT(float value, time_t timestamp);
+
+float MEOS_TFloatInst_value(MEOS_TFloatInst *tfloatinst);
+time_t MEOS_TFloatInst_timestamp(MEOS_TFloatInst *tfloatinst);
+char *MEOS_TFloatInst_str(MEOS_TFloatInst *tfloatinst);
+
+void MEOS_deleteTFloatInst(MEOS_TFloatInst *tfloatinst);
+
+/*****************************************************************************
+ * TFloatSeq
+ *****************************************************************************/
+
+// Constructor with a string
+MEOS_TFloatSeq *MEOS_newTFloatSeq(char *serialized);
+// Constructor with a set of instants, two boolean bounds and interpolation
+MEOS_TFloatSeq *MEOS_newTFloatSeq_IBBI(MEOS_TFloatInst **instants, int count, bool lower_inc,
+                                       bool upper_inc, MEOS_Interpolation interpolation);
+// Constructor with a set of instants as strings, two boolean bounds and interpolation
+MEOS_TFloatSeq *MEOS_newTFloatSeq_IsBBI(char **instants, int count, bool lower_inc, bool upper_inc,
+                                        MEOS_Interpolation interpolation);
+
+bool MEOS_TFloatSeq_lower_inc(MEOS_TFloatSeq *tfloatseq);
+bool MEOS_TFloatSeq_upper_inc(MEOS_TFloatSeq *tfloatseq);
+MEOS_Interpolation MEOS_TFloatSeq_interpolation(MEOS_TFloatSeq *tfloatseq);
+MEOS_TFloatInst **MEOS_TFloatSeq_instants(MEOS_TFloatSeq *tfloatseq, int *count);
+char *MEOS_TFloatSeq_str(MEOS_TFloatSeq *tfloatseq);
+
+void MEOS_deleteTFloatSeq(MEOS_TFloatSeq *tfloatseq);
+
+/*****************************************************************************
+ * TFloatSeqSet
+ *****************************************************************************/
+
+// Constructor with a string
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet(char *serialized);
+// Constructor with a set of sequences and interpolation
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet_SI(MEOS_TFloatSeq **sequences, int count,
+                                           MEOS_Interpolation interpolation);
+// Constructor with a set of sequences as strings and interpolation
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet_SsI(char **sequences, int count,
+                                            MEOS_Interpolation interpolation);
+
+MEOS_TFloatSeq **MEOS_TFloatSeqSet_sequences(MEOS_TFloatSeqSet *tfloatseqset, int *count);
+MEOS_TFloatSeqSet *MEOS_TFloatSeqSet_atPeriod(MEOS_TFloatSeqSet *tfloatseqset, MEOS_Period *period);
+char *MEOS_TFloatSeqSet_str(MEOS_TFloatSeqSet *tfloatseqset);
+
+void MEOS_deleteTFloatSeqSet(MEOS_TFloatSeqSet *tfloatseqset);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/capi/source/meos_c.cpp b/capi/source/meos_c.cpp
new file mode 100644
index 0000000..a3a4fe3
--- /dev/null
+++ b/capi/source/meos_c.cpp
@@ -0,0 +1,268 @@
+#include <meos/meos_c.h>
+
+#include <cstring>
+#include <ctime>
+#include <meos/types/box/TBox.hpp>
+#include <meos/types/temporal/TSequence.hpp>
+#include <meos/types/temporal/TSequenceSet.hpp>
+#include <meos/types/time/Period.hpp>
+#include <sstream>
+
+extern "C" {
+
+/*****************************************************************************
+ * TBox
+ *****************************************************************************/
+
+MEOS_TBox *MEOS_newTBox(double const xmin, double const xmax) {
+  return reinterpret_cast<MEOS_TBox *>(new meos::TBox(xmin, xmax));
+}
+
+double MEOS_TBox_xmin(MEOS_TBox *tbox) {
+  auto t = reinterpret_cast<meos::TBox *>(tbox);
+  return t->xmin();
+}
+
+double MEOS_TBox_xmax(MEOS_TBox *tbox) {
+  auto t = reinterpret_cast<meos::TBox *>(tbox);
+  return t->xmax();
+}
+
+void MEOS_deleteTBox(MEOS_TBox *tbox) {
+  auto t = reinterpret_cast<meos::TBox *>(tbox);
+  delete t;
+}
+
+/*****************************************************************************
+ * Period
+ *****************************************************************************/
+
+MEOS_Period *MEOS_newPeriod(char *serialized) {
+  return reinterpret_cast<MEOS_Period *>(new meos::Period(serialized));
+}
+
+MEOS_Period *MEOS_newPeriod_TT(time_t lower, time_t upper) {
+  return reinterpret_cast<MEOS_Period *>(
+      new meos::Period(std::chrono::system_clock::from_time_t(lower),
+                       std::chrono::system_clock::from_time_t(upper)));
+}
+
+MEOS_Period *MEOS_newPeriod_TTBB(time_t lower, time_t upper, bool lower_inc, bool upper_inc) {
+  return reinterpret_cast<MEOS_Period *>(
+      new meos::Period(std::chrono::system_clock::from_time_t(lower),
+                       std::chrono::system_clock::from_time_t(upper), lower_inc, upper_inc));
+}
+
+// Remember to free the result!
+char *MEOS_Period_str(MEOS_Period *period) {
+  auto t = reinterpret_cast<meos::Period *>(period);
+  std::stringstream output;
+  output << *t;
+  auto s = output.str();
+  const std::size_t len = s.length();
+  char *result = new char[len + 1];
+  if (result) {
+    std::strcpy(result, s.c_str());
+  }
+  return result;
+}
+
+void MEOS_deletePeriod(MEOS_Period *period) {
+  auto t = reinterpret_cast<meos::Period *>(period);
+  delete t;
+}
+
+/*****************************************************************************
+ * TFloatInst
+ *****************************************************************************/
+
+MEOS_TFloatInst *MEOS_newTFloatInst(char *serialized) {
+  return reinterpret_cast<MEOS_TFloatInst *>(new meos::TFloatInst(serialized));
+}
+
+MEOS_TFloatInst *MEOS_newTFloatInst_VT(float value, time_t timestamp) {
+  return reinterpret_cast<MEOS_TFloatInst *>(
+      new meos::TFloatInst(value, std::chrono::system_clock::from_time_t(timestamp)));
+}
+
+float MEOS_TFloatInst_value(MEOS_TFloatInst *tfloatinst) {
+  auto t = reinterpret_cast<meos::TFloatInst *>(tfloatinst);
+  return t->getValue();
+}
+
+time_t MEOS_TFloatInst_timestamp(MEOS_TFloatInst *tfloatinst) {
+  auto t = reinterpret_cast<meos::TFloatInst *>(tfloatinst);
+  return std::chrono::system_clock::to_time_t(t->getTimestamp());
+}
+
+// Remember to free the result!
+char *MEOS_TFloatInst_str(MEOS_TFloatInst *tfloatinst) {
+  auto t = reinterpret_cast<meos::TFloatInst *>(tfloatinst);
+  std::stringstream output;
+  output << *t;
+  auto s = output.str();
+  const std::size_t len = s.length();
+  char *result = new char[len + 1];
+  if (result) {
+    std::strcpy(result, s.c_str());
+  }
+  return result;
+}
+
+void MEOS_deleteTFloatInst(MEOS_TFloatInst *tfloatinst) {
+  auto t = reinterpret_cast<meos::TFloatInst *>(tfloatinst);
+  delete t;
+}
+
+/*****************************************************************************
+ * TFloatSeq
+ *****************************************************************************/
+
+MEOS_TFloatSeq *MEOS_newTFloatSeq(char *serialized) {
+  return reinterpret_cast<MEOS_TFloatSeq *>(new meos::TFloatSeq(serialized));
+}
+
+MEOS_TFloatSeq *MEOS_newTFloatSeq_IBBI(MEOS_TFloatInst **instants, int count, bool lower_inc,
+                                       bool upper_inc, MEOS_Interpolation interpolation) {
+  meos::TFloatInst *insts = new meos::TFloatInst[count];
+  for (size_t i = 0; i < count; i++) {
+    insts[i] = *reinterpret_cast<meos::TFloatInst *>(instants[i]);
+  }
+  std::set<meos::TFloatInst> s(insts, insts + count);
+  auto interp = interpolation == MEOS_Interpolation_Linear ? meos::Interpolation::Linear
+                                                           : meos::Interpolation::Stepwise;
+  return reinterpret_cast<MEOS_TFloatSeq *>(new meos::TFloatSeq(s, lower_inc, upper_inc, interp));
+}
+
+MEOS_TFloatSeq *MEOS_newTFloatSeq_IsBBI(char **instants, int count, bool lower_inc, bool upper_inc,
+                                        MEOS_Interpolation interpolation) {
+  std::set<meos::TFloatInst> s;
+  for (size_t i = 0; i < count; i++) {
+    meos::TFloatInst inst(instants[i]);
+    s.insert(inst);
+  }
+  auto interp = interpolation == MEOS_Interpolation_Linear ? meos::Interpolation::Linear
+                                                           : meos::Interpolation::Stepwise;
+  return reinterpret_cast<MEOS_TFloatSeq *>(new meos::TFloatSeq(s, lower_inc, upper_inc, interp));
+}
+
+bool MEOS_TFloatSeq_lower_inc(MEOS_TFloatSeq *tfloatseq) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  return t->lower_inc();
+}
+
+bool MEOS_TFloatSeq_upper_inc(MEOS_TFloatSeq *tfloatseq) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  return t->upper_inc();
+}
+
+MEOS_Interpolation MEOS_TFloatSeq_interpolation(MEOS_TFloatSeq *tfloatseq) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  return t->interpolation() == meos::Interpolation::Linear ? MEOS_Interpolation_Linear
+                                                           : MEOS_Interpolation_Stepwise;
+}
+
+MEOS_TFloatInst **MEOS_TFloatSeq_instants(MEOS_TFloatSeq *tfloatseq, int *count) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  auto r = t->instants();
+  *count = r.size();
+  MEOS_TFloatInst **instants = new MEOS_TFloatInst *[*count];
+  size_t i = 0;
+  for (auto const &it : r) {
+    auto x = new meos::TFloatInst(it);
+    instants[i++] = reinterpret_cast<MEOS_TFloatInst *>(x);
+  }
+  return instants;
+}
+
+// Remember to free the result!
+char *MEOS_TFloatSeq_str(MEOS_TFloatSeq *tfloatseq) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  std::stringstream output;
+  output << *t;
+  auto s = output.str();
+  const std::size_t len = s.length();
+  char *result = new char[len + 1];
+  if (result) {
+    std::strcpy(result, s.c_str());
+  }
+  return result;
+}
+
+void MEOS_deleteTFloatSeq(MEOS_TFloatSeq *tfloatseq) {
+  auto t = reinterpret_cast<meos::TFloatSeq *>(tfloatseq);
+  delete t;
+}
+
+/*****************************************************************************
+ * TFloatSeqSet
+ *****************************************************************************/
+
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet(char *serialized) {
+  return reinterpret_cast<MEOS_TFloatSeqSet *>(new meos::TFloatSeqSet(serialized));
+}
+
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet_SI(MEOS_TFloatSeq **sequences, int count,
+                                           MEOS_Interpolation interpolation) {
+  meos::TFloatSeq *seqs = new meos::TFloatSeq[count];
+  for (size_t i = 0; i < count; i++) {
+    seqs[i] = *reinterpret_cast<meos::TFloatSeq *>(sequences[i]);
+  }
+  std::set<meos::TFloatSeq> s(seqs, seqs + count);
+  auto interp = interpolation == MEOS_Interpolation_Linear ? meos::Interpolation::Linear
+                                                           : meos::Interpolation::Stepwise;
+  return reinterpret_cast<MEOS_TFloatSeqSet *>(new meos::TFloatSeqSet(s, interp));
+}
+
+MEOS_TFloatSeqSet *MEOS_newTFloatSeqSet_SsI(char **sequences, int count,
+                                            MEOS_Interpolation interpolation) {
+  std::set<meos::TFloatSeq> s;
+  for (size_t i = 0; i < count; i++) {
+    meos::TFloatSeq seq(sequences[i]);
+    s.insert(seq);
+  }
+  auto interp = interpolation == MEOS_Interpolation_Linear ? meos::Interpolation::Linear
+                                                           : meos::Interpolation::Stepwise;
+  return reinterpret_cast<MEOS_TFloatSeqSet *>(new meos::TFloatSeqSet(s, interp));
+}
+
+MEOS_TFloatSeq **MEOS_TFloatSeqSet_sequences(MEOS_TFloatSeqSet *tfloatseqset, int *count) {
+  auto t = reinterpret_cast<meos::TFloatSeqSet *>(tfloatseqset);
+  auto r = t->sequences();
+  *count = r.size();
+  MEOS_TFloatSeq **sequences = new MEOS_TFloatSeq *[*count];
+  size_t i = 0;
+  for (auto const &it : r) {
+    auto x = new meos::TFloatSeq(it);
+    sequences[i++] = reinterpret_cast<MEOS_TFloatSeq *>(x);
+  }
+  return sequences;
+}
+
+MEOS_TFloatSeqSet *MEOS_TFloatSeqSet_atPeriod(MEOS_TFloatSeqSet *tfloatseqset,
+                                              MEOS_Period *period) {
+  auto t = reinterpret_cast<meos::TFloatSeqSet *>(tfloatseqset);
+  auto p = reinterpret_cast<meos::Period *>(period);
+  auto r = t->atPeriod(*p);
+  return reinterpret_cast<MEOS_TFloatSeqSet *>(new meos::TFloatSeqSet(r));
+}
+
+// Remember to free the result!
+char *MEOS_TFloatSeqSet_str(MEOS_TFloatSeqSet *tfloatseqset) {
+  auto t = reinterpret_cast<meos::TFloatSeqSet *>(tfloatseqset);
+  std::stringstream output;
+  output << *t;
+  auto s = output.str();
+  const std::size_t len = s.length();
+  char *result = new char[len + 1];
+  if (result) {
+    std::strcpy(result, s.c_str());
+  }
+  return result;
+}
+
+void MEOS_deleteTFloatSeqSet(MEOS_TFloatSeqSet *tfloatseqset) {
+  auto t = reinterpret_cast<meos::TFloatSeqSet *>(tfloatseqset);
+  delete t;
+}
+}
diff --git a/include/meos/types/temporal/TSequence.hpp b/include/meos/types/temporal/TSequence.hpp
index 91524a4..13c3b62 100644
--- a/include/meos/types/temporal/TSequence.hpp
+++ b/include/meos/types/temporal/TSequence.hpp
@@ -62,6 +62,7 @@ template <typename BaseType = float> class TSequence : public TemporalSet<BaseTy
   TSequence<BaseType> *shift_impl(duration_ms const timedelta) const override;
   bool intersectsTimestamp(time_point const datetime) const override;
   bool intersectsPeriod(Period const period) const override;
+  TSequence<BaseType> atPeriod(Period const &period) const;
 
   std::istream &read(std::istream &in, bool with_interp = true, bool with_srid = true);
   std::ostream &write(std::ostream &os, bool with_interp = true, bool with_srid = true) const;
@@ -110,6 +111,9 @@ template <typename BaseType = float> class TSequence : public TemporalSet<BaseTy
   std::ostream &write_internal(std::ostream &os, bool with_interp = true) const;
 
   TSequence<BaseType> *clone_impl() const override { return new TSequence<BaseType>(*this); };
+
+  TInstant<BaseType> atTimestamp(time_point t) const;
+  int findTimestamp(time_point t) const;
 };
 
 typedef TSequence<bool> TBoolSeq;
diff --git a/include/meos/types/temporal/TSequenceSet.hpp b/include/meos/types/temporal/TSequenceSet.hpp
index fd9bf20..4e17e5c 100644
--- a/include/meos/types/temporal/TSequenceSet.hpp
+++ b/include/meos/types/temporal/TSequenceSet.hpp
@@ -93,6 +93,7 @@ template <typename BaseType = float> class TSequenceSet
   TSequenceSet<BaseType> *shift_impl(duration_ms const timedelta) const override;
   bool intersectsTimestamp(time_point const datetime) const override;
   bool intersectsPeriod(Period const period) const override;
+  TSequenceSet<BaseType> atPeriod(Period const &period) const;
 
   std::istream &read(std::istream &in);
   std::ostream &write(std::ostream &os) const;
@@ -137,6 +138,8 @@ template <typename BaseType = float> class TSequenceSet
   std::ostream &write_internal(std::ostream &os) const;
 
   TSequenceSet<BaseType> *clone_impl() const override { return new TSequenceSet<BaseType>(*this); };
+
+  int findTimestamp(time_point t) const;
 };
 
 typedef TSequenceSet<bool> TBoolSeqSet;
diff --git a/include/meos/types/temporal/Temporal.hpp b/include/meos/types/temporal/Temporal.hpp
index 445834f..5103a15 100644
--- a/include/meos/types/temporal/Temporal.hpp
+++ b/include/meos/types/temporal/Temporal.hpp
@@ -133,6 +133,17 @@ template <typename BaseType = float> class Temporal
    */
   bool intersectsPeriodSet(PeriodSet const periodset) const;
 
+  /**
+   * @brief Restrict to a period
+   * @code
+   * TFloatSeqSet t("{[1@2012-01-01, 3@2012-01-03), [3@2012-01-04, 1@2012-01-06)}");
+   * Period p("[2012-01-02,2012-01-05)");
+   *
+   * t.atPeriod(p);  // <TFloatSeqSet {[2@2012-01-02, 3@2012-01-03), [3@2012-01-04, 2@2012-01-05)}>
+   * @endcode
+   */
+  // virtual BaseType atPeriod(Period period) const = 0;
+
 private:
   virtual Temporal<BaseType> *clone_impl() const = 0;
 
diff --git a/include/meos/types/time/Period.hpp b/include/meos/types/time/Period.hpp
index f44c1fb..6829fbf 100644
--- a/include/meos/types/time/Period.hpp
+++ b/include/meos/types/time/Period.hpp
@@ -43,7 +43,9 @@ class Period {
   duration_ms timespan() const;
   std::unique_ptr<Period> shift(duration_ms const timedelta) const;
   bool overlap(Period const &period) const;
+  bool contains(Period const &period) const;
   bool contains_timestamp(time_point const timestamp) const;
+  Period intersection(Period const &period) const;
 
   friend bool operator==(Period const &lhs, Period const &rhs);
   friend bool operator!=(Period const &lhs, Period const &rhs);
diff --git a/source/types/temporal/TSequence.cpp b/source/types/temporal/TSequence.cpp
index c253128..25a0576 100644
--- a/source/types/temporal/TSequence.cpp
+++ b/source/types/temporal/TSequence.cpp
@@ -182,30 +182,31 @@ int TSequence<BaseType>::compare_internal(Temporal<BaseType> const &other) const
   }
 
   TSequence<BaseType> const *that = dynamic_cast<TSequence<BaseType> const *>(&other);
-  // Compare number of instants
-  if (this->m_instants.size() < that->m_instants.size()) return -1;
-  if (this->m_instants.size() > that->m_instants.size()) return 1;
-
-  // Compare bounds
-  // [ < (, ) < ]
-  if ((this->m_lower_inc && !that->m_lower_inc) || (!this->m_upper_inc && that->m_upper_inc))
-    return -1;
-  // ( > [, ] > )
-  if ((that->m_lower_inc && !this->m_lower_inc) || (!that->m_upper_inc && this->m_upper_inc))
-    return 1;
 
   // Compare instant by instant
   auto lhs_instants = this->instants();
   auto rhs_instants = that->instants();
   auto lhs = lhs_instants.begin();
   auto rhs = rhs_instants.begin();
-  while (lhs != lhs_instants.end()) {
+  while (lhs != lhs_instants.end() && rhs != rhs_instants.end()) {
     if (*lhs < *rhs) return -1;
     if (*lhs > *rhs) return 1;
     lhs++;
     rhs++;
   }
 
+  // Compare number of instants
+  if (this->m_instants.size() < that->m_instants.size()) return -1;
+  if (this->m_instants.size() > that->m_instants.size()) return 1;
+
+  // Compare bounds
+  // [ < (, ) < ]
+  if ((this->m_lower_inc && !that->m_lower_inc) || (!this->m_upper_inc && that->m_upper_inc))
+    return -1;
+  // ( > [, ] > )
+  if ((that->m_lower_inc && !this->m_lower_inc) || (!that->m_upper_inc && this->m_upper_inc))
+    return 1;
+
   // Compare Interpolation
   if (this->interpolation() < that->interpolation()) return -1;
   if (this->interpolation() > that->interpolation()) return 1;
@@ -298,6 +299,164 @@ template <typename BaseType> bool TSequence<BaseType>::intersectsPeriod(Period c
   return this->period().overlap(period);
 }
 
+template <typename BaseType>
+TSequence<BaseType> TSequence<BaseType>::atPeriod(Period const &period) const {
+  Period inter = this->period().intersection(period);
+  auto interp = this->interpolation();
+
+  /* Intersecting period is instantaneous */
+  if (inter.lower() == inter.upper()) {
+    set<TInstant<BaseType>> instants;
+    TInstant<BaseType> instant = this->atTimestamp(inter.lower());
+    instants.insert(instant);
+    return TSequence<BaseType>(instants, true, true, interp);
+  }
+
+  int n = this->findTimestamp(inter.lower());
+
+  /* If the lower bound of the intersecting period is exclusive */
+  if (n == -1) n = 0;
+  set<TInstant<BaseType>> instants;
+
+  /* Compute the value at the beginning of the intersecting period */
+  TInstant<BaseType> inst1 = this->instantN(n);
+  TInstant<BaseType> inst2 = this->instantN(n + 1);
+  instants.insert(at_timestamp1(inst1, inst2, interp, inter.lower()));
+  for (size_t i = n + 2; i < this->m_instants.size(); i++) {
+    /* If the end of the intersecting period is between inst1 and inst2 */
+    if (inst1.getTimestamp() <= inter.upper() && inter.upper() <= inst2.getTimestamp()) break;
+
+    inst1 = inst2;
+    inst2 = this->instantN(i);
+    /* If the intersecting period contains inst1 */
+    if (inter.lower() <= inst1.getTimestamp() && inst1.getTimestamp() <= inter.upper())
+      instants.insert(inst1);
+  }
+
+  /* The last two values of sequences with step interpolation and
+   * exclusive upper bound must be equal */
+  if (interp == Interpolation::Linear || inter.upper_inc())
+    instants.insert(at_timestamp1(inst1, inst2, interp, inter.upper()));
+  else {
+    auto value = instants.rbegin()->getValue();
+    instants.insert(TInstant<BaseType>(value, inter.upper()));
+  }
+
+  /* Since by definition the sequence is normalized it is not necessary to
+   * normalize the projection of the sequence to the period */
+  TSequence<BaseType> result(instants, inter.lower_inc(), inter.upper_inc(), interp);
+  return result;
+}
+
+template <typename BaseType> int TSequence<BaseType>::findTimestamp(time_point t) const {
+  int first = 0;
+  int const seq_size = this->m_instants.size();
+  int last = seq_size - 2;
+  int middle = (first + last) / 2;
+  while (first <= last) {
+    TInstant<BaseType> inst1 = this->instantN(middle);
+    TInstant<BaseType> inst2 = this->instantN(middle + 1);
+    bool lower_inc = (middle == 0) ? this->period().lower_inc() : true;
+    bool upper_inc = (middle == seq_size - 2) ? this->period().upper_inc() : false;
+    if ((inst1.getTimestamp() < t && t < inst2.getTimestamp())
+        || (lower_inc && inst1.getTimestamp() == t) || (upper_inc && inst2.getTimestamp() == t))
+      return middle;
+    if (t <= inst1.getTimestamp())
+      last = middle - 1;
+    else
+      first = middle + 1;
+    middle = (first + last) / 2;
+  }
+  return -1;
+}
+
+/**
+ * Returns a point interpolated from the geometry/geography segment with
+ * respect to the fraction of its total length.
+ *
+ * @param[in] start,end Points defining the segment
+ * @param[in] ratio Float between 0 and 1 representing the fraction of the
+ * total length of the segment where the point must be located
+ */
+GeomPoint geoseg_interpolate_point(GeomPoint start, GeomPoint end, double ratio) {
+  int srid = start.srid();
+  // TODO add support for geodetic
+  // bool geodetic = start.geodetic();
+  GeomPoint p(start.x() + ((end.x() - start.x()) * ratio),
+              start.y() + ((end.y() - start.y()) * ratio), srid);
+  return p;
+}
+
+template <typename BaseType> BaseType value_at_timestamp1(TInstant<BaseType> const & /*inst1*/,
+                                                          TInstant<BaseType> const & /*inst2*/,
+                                                          Interpolation /*interpolation*/,
+                                                          time_point /*t*/) {
+  // Check template specializations for double and GeomPoint
+  throw invalid_argument("Unsupported base type for function");
+}
+
+template <> float value_at_timestamp1(TInstant<float> const &inst1, TInstant<float> const &inst2,
+                                      Interpolation interpolation, time_point t) {
+  auto start = inst1.getValue();
+  auto end = inst2.getValue();
+  /* Constant segment or t is equal to lower bound or step interpolation */
+  if (start == end || inst1.getTimestamp() == t
+      || (interpolation != Interpolation::Linear && t < inst2.getTimestamp()))
+    return start;
+
+  /* t is equal to upper bound */
+  if (inst2.getTimestamp() == t) return end;
+
+  /* Interpolation for types with linear interpolation */
+  auto duration1 = t - inst1.getTimestamp();
+  auto duration2 = inst2.getTimestamp() - inst1.getTimestamp();
+  float ratio = ((float)duration1.count()) / duration2.count();
+  return start + (end - start) * ratio;
+}
+
+template <> GeomPoint value_at_timestamp1(TInstant<GeomPoint> const &inst1,
+                                          TInstant<GeomPoint> const &inst2,
+                                          Interpolation interpolation, time_point t) {
+  auto start = inst1.getValue();
+  auto end = inst2.getValue();
+  /* Constant segment or t is equal to lower bound or step interpolation */
+  if (start == end || inst1.getTimestamp() == t
+      || (interpolation != Interpolation::Linear && t < inst2.getTimestamp()))
+    return start;
+
+  /* t is equal to upper bound */
+  if (inst2.getTimestamp() == t) return end;
+
+  /* Interpolation for types with linear interpolation */
+  auto duration1 = t - inst1.getTimestamp();
+  auto duration2 = inst2.getTimestamp() - inst1.getTimestamp();
+  double ratio = duration1 / duration2;
+  return geoseg_interpolate_point(start, end, ratio);
+}
+
+template <typename BaseType>
+TInstant<BaseType> at_timestamp1(TInstant<BaseType> const &inst1, TInstant<BaseType> const &inst2,
+                                 Interpolation interpolation, time_point t) {
+  auto value = value_at_timestamp1(inst1, inst2, interpolation, t);
+  return TInstant<BaseType>(value, t);
+}
+
+template <typename BaseType>
+TInstant<BaseType> TSequence<BaseType>::atTimestamp(time_point t) const {
+  /* Bounding box test */
+  if (!this->period().contains_timestamp(t))
+    throw invalid_argument("Specified timestamp does not overlap with this sequence");
+
+  /* Instantaneous sequence */
+  if (this->m_instants.size() == 1) return this->startInstant();
+
+  /* General case */
+  int n = this->findTimestamp(t);
+  TInstant<BaseType> inst1 = this->instantN(n);
+  TInstant<BaseType> inst2 = this->instantN(n + 1);
+  return at_timestamp1(inst1, inst2, this->interpolation(), t);
+}
+
 template <typename BaseType>
 istream &TSequence<BaseType>::read_internal(istream &in, bool with_interp) {
   char c;
diff --git a/source/types/temporal/TSequenceSet.cpp b/source/types/temporal/TSequenceSet.cpp
index fd806d9..6aa92c1 100644
--- a/source/types/temporal/TSequenceSet.cpp
+++ b/source/types/temporal/TSequenceSet.cpp
@@ -337,6 +337,53 @@ bool TSequenceSet<BaseType>::intersectsPeriod(Period const period) const {
   return false;
 }
 
+template <typename BaseType>
+TSequenceSet<BaseType> TSequenceSet<BaseType>::atPeriod(Period const &period) const {
+  int loc = this->findTimestamp(period.lower());
+  int const seq_size = this->m_sequences.size();
+  if (loc >= seq_size) {
+    throw invalid_argument("no value at specified period");
+  }
+  TSequence<BaseType> seq;
+  set<TSequence<BaseType>> sequences;
+  for (int i = loc; i < seq_size; i++) {
+    seq = this->sequenceN(i);
+    auto seq_period = seq.period();
+    if (period.contains(seq_period)) {
+      sequences.insert(seq);
+    } else if (period.overlap(seq_period)) {
+      sequences.insert(seq.atPeriod(period));
+    }
+    if ((period.upper() < seq_period.upper())
+        || ((period.upper() == seq_period.upper()) && seq_period.upper_inc())) {
+      break;
+    }
+  }
+  if (sequences.size() == 0) {
+    throw invalid_argument("no value at specified period");
+  }
+  return TSequenceSet<BaseType>(sequences);
+}
+
+template <typename BaseType> int TSequenceSet<BaseType>::findTimestamp(time_point t) const {
+  int first = 0, middle = 0, last = this->m_sequences.size() - 1;
+  TSequence<BaseType> seq;
+  while (first <= last) {
+    middle = (first + last) / 2;
+    seq = this->sequenceN(middle);
+    Period p = seq.period();
+    if (p.contains_timestamp(t)) {
+      return middle;
+    }
+    if (t <= p.lower())
+      last = middle - 1;
+    else
+      first = middle + 1;
+  }
+  if (t >= seq.period().upper()) middle++;
+  return middle;
+}
+
 template <typename BaseType> istream &TSequenceSet<BaseType>::read_internal(istream &in) {
   char c;
 
diff --git a/source/types/time/Period.cpp b/source/types/time/Period.cpp
index 5d0601c..f4c996c 100644
--- a/source/types/time/Period.cpp
+++ b/source/types/time/Period.cpp
@@ -84,12 +84,38 @@ bool Period::overlap(Period const &period) const {
                                         : period.upper_inc() && this->lower_inc();
 }
 
+bool Period::contains(Period const &period) const {
+  bool c1 = this->lower() < period.lower();
+  bool c2 = this->upper() > period.upper();
+  if ((c1 || (this->lower() == period.lower() && (this->lower_inc() || !period.lower_inc())))
+      && (c2 || (this->upper() == period.upper() && (this->upper_inc() || !period.upper_inc())))) {
+    return true;
+  }
+  return false;
+}
+
 bool Period::contains_timestamp(time_point const timestamp) const {
   return ((this->lower() < timestamp && timestamp < this->upper())
           || (this->lower_inc() && this->lower() == timestamp)
           || (this->upper_inc() && this->upper() == timestamp));
 }
 
+Period Period::intersection(Period const &other) const {
+  /* Bounding box test */
+  if (!this->overlap(other))
+    throw invalid_argument("Specified period does not overlap with this period");
+
+  time_point lower = max(this->lower(), other.lower());
+  time_point upper = min(this->upper(), other.upper());
+  bool lower_inc = this->lower() == other.lower()
+                       ? this->lower_inc() && other.lower_inc()
+                       : (lower == this->lower() ? this->lower_inc() : other.lower_inc());
+  bool upper_inc = this->upper() == other.upper()
+                       ? this->upper_inc() && other.upper_inc()
+                       : (upper == this->upper() ? this->upper_inc() : other.upper_inc());
+  return Period(lower, upper, lower_inc, upper_inc);
+}
+
 int Period::compare(Period const &other) const {
   if (lower() < other.lower())
     return -1;
diff --git a/test/source/types/temporal/tsequence.cpp b/test/source/types/temporal/tsequence.cpp
index e6d7719..e2a717e 100644
--- a/test/source/types/temporal/tsequence.cpp
+++ b/test/source/types/temporal/tsequence.cpp
@@ -216,8 +216,8 @@ TEMPLATE_TEST_CASE("TSequence comparision operators", "[tsequence]", int, float)
   SECTION("lhs < rhs") {
     SECTION("different sizes") {
       set<TInstant<TestType>> lhs_instants = {
-          TInstant<TestType>(4, unix_time_point(2012, 1, 4)),
-          TInstant<TestType>(5, unix_time_point(2012, 1, 5)),
+          TInstant<TestType>(1, unix_time_point(2012, 1, 1)),
+          TInstant<TestType>(2, unix_time_point(2012, 1, 2)),
       };
       TSequence<TestType> lhs(lhs_instants, true, false);
       set<TInstant<TestType>> rhs_instants = {
diff --git a/test/source/types/temporal/tsequenceset.cpp b/test/source/types/temporal/tsequenceset.cpp
index 2d43487..438b1a4 100644
--- a/test/source/types/temporal/tsequenceset.cpp
+++ b/test/source/types/temporal/tsequenceset.cpp
@@ -576,3 +576,25 @@ TEST_CASE("interpolation is maintained in sequences", "[tsequenceset]") {
     REQUIRE(seq.interpolation() == Interpolation::Stepwise);
   }
 }
+
+TEST_CASE("atPeriod", "[tsequenceset]") {
+  SECTION("no overlap at sequence level") {
+    TFloatSeqSet seqset("{[1@2012-01-02, 3@2012-01-03), [3@2012-01-04, 1@2012-01-06)}");
+    Period period("[2012-01-01,2012-01-07)");
+    TFloatSeqSet result = seqset.atPeriod(period);
+
+    stringstream output;
+    output << result;
+    REQUIRE(output.str() == "{[1@2012-01-02T00:00:00+0000, 3@2012-01-03T00:00:00+0000), [3@2012-01-04T00:00:00+0000, 1@2012-01-06T00:00:00+0000)}");
+  }
+
+  SECTION("overlap at sequence level") {
+    TFloatSeqSet seqset("{[1@2012-01-01, 3@2012-01-03), [3@2012-01-04, 1@2012-01-06)}");
+    Period period("[2012-01-02,2012-01-05)");
+    TFloatSeqSet result = seqset.atPeriod(period);
+
+    stringstream output;
+    output << result;
+    REQUIRE(output.str() == "{[2@2012-01-02T00:00:00+0000, 3@2012-01-03T00:00:00+0000), [3@2012-01-04T00:00:00+0000, 2@2012-01-05T00:00:00+0000)}");
+  }
+}