The gradual decommissioning of fossil fuel-driven power plants, that traditionally provide most operational flexibility in power systems, has led to more frequent grid stability issues. To compensate for the lack of flexible resources, Distributed Energy Resources (DERs) in distribution networks can be employed. To facilitate the use of DERs, the aggregated flexibility in a distribution grid is commonly represented on a PQ-plane displaying the feasible active and reactive power exchange with the upstream grid. This paper proposes a fast feasible operating region mapping mechanism that utilizes a linear power flow approximation in combination with linearized generator, current, and voltage constraints to construct a high-dimensional polyhedral feasible set of DER power injections. The obtained polytope is projected onto the PQ-plane using Fourier-Motzkin Elimination to represent the aggregate network flexibility. Additionally, uncertainty in DER generation is addressed using chance-constraints. Our analysis of a modified 33-bus IEEE test system demonstrates that the proposed method obtains more accurate approximations than former geometric methods and is ten times faster than the tested optimization-based method.