Many scientific phenomena are modeled by Partial Differential Equations (PDEs). The development of data gathering tools along with the advances in machine learning (ML) techniques have raised opportunities for data-driven identification of governing equations from experimentally observed data. We propose an ML method to discover the terms involved in the equation from two-dimensional spatiotemporal data. Robust and useful physical features are extracted from data samples to represent the specific behaviors imposed by each mathematical term in the equation. Compared to the previous models, this idea provides us with the ability to discover 2D equations with time derivatives of different orders, and also to identify new underlying physics on which the model has not been trained. Moreover, the model can work with small sets of low-resolution data while avoiding numerical differentiations. The results indicate robustness of the features extracted based on prior knowledge in comparison to automatically detected features by a Three-dimensional Convolutional Neural Network (3D CNN) given the same amounts of data. Although particular PDEs are studied in this work, the idea of the proposed approach could be extended for reliable identification of various PDEs.