The Martian isochrons are the basis in analyzing the impact flux and surface dating, and currently, they are usually derived from those of the Moon because no samples have been collected from Mars. However, the isochrons obtained by this method have substantial uncertainty, and they urgently need to be optimized based on the samples that are about to be obtained. To support the upcoming Mars sample return missions, we utilized a high-resolution Context Camera mosaic and selected 17 regions from diverse geological units and ages across Mars's surface to establish an observed Martian crater production function (PF). Craters were manually mapped across these regions. The crater size-frequency distributions (CSFDs) from these regions have a strong correlation at the same diameter range on a logarithmic scale, indicating that they share a similar distribution shape regardless of geological units and ages. We obtained 155 effective CSFD bins suitable for fitting the crater PF. After testing on different fitting functions, we finally obtained the crater PF for the Martian surface over the diameter range of 0.15 to 13.5 km. There were significant differences between the directly mapped Mars crater PF and those derived from lunar models. In addition, the CSFDs obtained by previous researchers when doing dating works on the Martian surface are more consistent with the newly established crater PF than with the earlier proposed PFs. With the radiometric ages of the samples returned by future Mars sample return missions, this research could become the basis for establishing a new chronology system for Mars.