In the late 1980s, the lattice-Boltzmann method (LBM) emerged as an alternative approach to computational fluid dynamics. Instead of discretizing the Navier-Stokes equations directly, LBM discretizes the underlying Boltzmann equation. This leads to a local, explicit formulation on a Cartesian lattice, making it well-suited for massive parallelization. However, complex geometries create arbitrarily cut cells in the Cartesian lattice, complicating boundary condition enforcement. We seek a boundary treatment that is stable, second-order accurate and conservative. The current state-of-the-art boundary treatment (volumetric) is stable and conservative, but not second order accurate. They take a finite-volume approach, assuming constant concentrations in each cell when evaluating the interaction with the boundary. We propose to further increase accuracy and promote stability by locally constructing piecewise polynomials of the populations. Boundary interactions are evaluated on these, the result is integrated back to cell lattice values. A 2D Python implementation of this approach is developed. It is first tested on a channel which is inclined relative to the lattice, producing arbitrarily cut cells. The velocity error shows second-order convergence. Furthermore, stability is verified through a novel eigenvalue analysis of the global streaming operation. Next, flow around a NACA0012 airfoil is tested and compared with a direct numerical simulation of the Navier-Stokes equations. We observe good agreement in the total forces. However, in the surface forces, some oscillations are observed with coarse meshes, which are common in methods based on a Cartesian grids. In the future, we intend to address the surface force oscillations, incorporate wall models, extend the implementation to 3D and parallelize it such that it becomes applicable to real-world cases.