Multi-criteria GIS analysis of the topography of the Moon and better solutions for potential landing

During the past twenty years, the need to reach the Moon by the private spacemissions has been growing. Some of the private missions are supported by Google Lunar X-prize and Space-X. In the period between 2020 and 2050 private companies will be planning landing to the Moon with their own capacity. These missions can send new geodesy and cartography data. Lunar topography modelling with new satellite and remote sensing data gives plenty of possibilities for its exploration. GIS (Geographical Information System) may be successfully to the Moon topography analysis. According to the results after GIS numerical analysis 30% of the territory of the Moon showed excellent characteristics for landing. The most useful parts of the Moon for potential landing belong to the altitude between 2,000 and 3,000 m and on the plateaus with the north-east direction of azimuth. These plateaus have an excellent inclination of 3∘ and azimuth of 120∘. The main aim and goal of this investigation would be in better understanding of Moon topography and relief. With help of GIS numerical methods, the astronomical geodesy may be applied in better way. A potential mission to the Moon can use this topography investigation, presented maps and results.


Introduction
In 1969, humankind reached the Moon. It had been a small step for the man but a giant leap for the mankind. After that mission, we had exactly six missions more, as NASA successfully launched seven Apollo missions to the Moon (Apollo 11-17). After 1972 we haven't had any human mission to the Moon or soft landing on the Moon surface. In fact, we have a pause for more than 35 years. However, humankind in an attempt to derive the most precise model of the Moon. Some papers also included Moon topography assessment. Michael et al. (2018), investigated early topography of Moon and bombarding of meteors in the geological past of the Moon. Many models of Moon topography can be applied to find ice influence on the surface. These analyses always included precise geodesy and cartography methods (Rubanenko and Aharonson 2017). With help of special ray system of mapping it is possible to find relationship into the primary and secondary (or measured craters), Elliott et al. (2018). The Moon is a body that orbits the planet Earth like its natural satellite. The Moon is the fifth natural satellite in the whole Solar system, when it comes to its volume.
In one important research about craters impacts the main goal was in determining between the craters caused by asteroid and comet impacts, the aperiodic component could be removed. The lunar impact craters may eventually provide the required supplementary data or new investigation (Jetsu 2017). The Moon was created before 4.51×10 9 years. Some scientists think that the Moon belongs to the Earth like one space body with co-formation of the Earth and the Moon together in the primordial accretion disk. The Moon's spin and orbital periods became locked together very early in its existence, when it was much closer to the Earth than it is now. The far side is, in average, about 5 km higher, with respect to the Moon's center of mass, then the near side. This hypothesis is not fully accepted. The far side of the Moon was formed by crust with the depth of layer of 48 km. The near side of the Moon was formed in a very similar way as the far side, but with deeper layer of 76 km Wilhelms et al. (1987); Fielder (1961); Ferrari et al. (1978); Binzel and Van Flandern (1979);De Hon and Waskom (1976). Since the near side is lower, volcanic magma has more easily found its way to the surface here, pouring from volcanic fissures into the low regions of the largest craters and solidifying to form the lunar seas. Surface geology of the Moon or topography of the Moon, in the last decades, has been measured with laser altimetry and techniques of image analysis. The most dominant terrain-geomorphologic structure is South Pole Aitken basin with the diameter of 2,240 km, which is also the largest crater and the lowest point on the Moon. The massive mountains' systems on the Moon are located on the north-east. The dark is called Maria (Latin Abbrev.) in transcription Sea. These basins are filled with lava and basalts. 31% on the near side of the Moon belong to the Sea area. Most of the Moon's mare basalts erupted during 3.0 to 3.6 billion years ago Ershkov (2015); Olson and Wilhems (1974); Matzger and Parker (1974). The second terrain objects are the impact craters, which were created by the impact of comets and asteroids. There are, according to the estimate, more than 300,000 craters on the near side of the Moon. These craters are wider than 1 km of diameter (see Figure 1). The Moon's geological timescale can be divided into three most prominent events involving creation of Nectaris, Imbrium and Orientale structures. Since the Moon has no atmosphere, we can perceive geological processes marked on the craters. Apollo missions gave results of the main stratigraphic structure between 3.8 to 4.1 billion years old. The thickness of the first layer called regolith is approximately 10 to 120 cm deep. Researchers in one investigation are established anaglyphonosphere axiomatic and landscape concepts obtained with the use of the axiomatic method. The first concept depicts the geographic envelope of the Moon as an anaglyphonosphere layer (relief) that is a continuum (total environment). (Kyryliuk and Kholiavchuk 2017).

LOLA (Topography, Surface Roughness, Derived Slope)
LOLA is the Lunar Orbiter Laser Altimeter, an instrument onboard LRO, that determines the global topography of the lunar surface at high resolution. LOLA data sets are produced by the LOLA Science Team at NASA Goddard Space Flight Center in Greenbelt, MD and numerous locations across the US. LOLA produces topographic data by firing a laser at the lunar surface and detecting the two-way travel time for the signal. LOLA's single laser, which is fired 28 times per second, is split into five beams by a diffractive optical element. LOLA sends a total of 140 pulses to the lunar surface each second, which enables the instrument to gather high resolution topographic data. This particular dataset shows the Moon's topography. The amount of time it takes for the laser pulses to return to LOLA's receiver onboard the LRO spacecraft indicates how far it is from the lunar surface, which can be used to create a topographic map of the lunar surface. The LOLA instrument onboard LRO has measured the altitude of more than 6 billion points on the lunar surface. Gray colors indicate high altitudes and blue and white indicate low altitudes ( Figure 1).

Methods and data
The Lunar Orbiter Laser Altimeter (LOLA) is producing a high-resolution global topographic model and geodetic framework that enables precise targeting, safe landing, and surface mobility to carry out exploratory activities. LOLA characterizes the polar illumination environment, and im-   Chapman et al. 1979;Eberhardt et al. 1970). The DEM and shape map of the (radius) of the Moon at a resolution of 128 pixels per degree, based on altimetry data acquired in September 2011 by the LOLA instrument, were used for GIS calculations. The ground tracks were interpolated using the Generic Mapping Tools programs -"surface" and "grid blend". The map is in the form of a binary table, with one row for every 0.007813 degrees of latitude, pixel registered. The map has covered the whole surface of the Moon within the radius of 1737400 m. Hill shade relief is enabled to give results in the GIS calculations. We used calculations and specially created algorithm to find, with the help of GIS multi-criteria decision analysis, potential places on the Moon for safe landing.

The main features on the Moon
There are two types of terrain on the Moon. The first type are dark plains, called maria (Latin for "sea"), the second terrain is cratered highland terrain. The whole surface was initially covered with craters during the time of massive bombardment in the past. Lava rose to the surface through cracks and fissures, filling the lower parts of the large craters to produce the dark plains. The plains reflect only 4% of the sunlight, whereas mountains reflect 11%. (Nash and Conel 1973).

Aristarchus Crater
This crater belongs to impact craters and it was created 300 million years ago with the diameter of 37 km. This young crater has a series of nested terraces which were produced by concentric slices of rock in the wall.

Mare (Sea) Crisium
Mare Crisium belongs to lava-filed impact crater sea, 3.9 billion years old. The diameter of this sea is 563 km. Mare Crisium has an extremely smooth floor, which varies in height by less than 90 m. Mare Crisium can be seen from the Earth, because over 95% of lunar craters are compactly circular.

Montes Apenninus (Mounts Apennines)
This object belongs to the mountain range and it is 3.9 billion years old. Length of this mountain is 401 km. The lunar Apennines form a ring around the southwestern edge of the Mare Imbrium impact basin.

Copernicus Crater
This crater belongs to an impact crater aged 900 million years and with the diameter of 91 km. This young ray crater has massive terraced walls.

Rupes Altai
Type of this object is a cliff, aged approximately 4.2 billion years. Its length is 507 km. Altai is, by far, the longest cliff on the Moon and about 1.8 km high, created due to the massive energy of the impact.

Alphonsus Crater
This object is of impact craters' type. It is 4 billion years old and 117 km long. The crater was formed in an impact, with the dark patches and fractures all around its floor.

Humboldt Crater
This crater also belongs to the impact type of craters. It is 3.8 billion years old and 189 km long. This crater is remarkable because its lava-filled floor is crisscrossed with a series of radial and concentric fractures.

Tycho Crater
This crater belongs to the impact type of craters. It is 100 million years old with the diameter of 85 km. It is lying in the southern highlands. Tycho is one of the most perfect, walled craters on the Moon, with a central mountain peak towering 3 km above a rough, infilled inner region.

GIS numerical criteria analysis
GIS and geo-data modelling in combination with new methods such as shading and subpixel analysis can present very interesting and powerful tool for image processing. The time of pixel and GIS numerical analysis started at the end of nineties. With the development of techniques and remote sensing, science' need for image processing is growing. Today, remote sensing and GIS analysis are not applied only to the lands, but they can as well be very successfully applied to the other celestial objects. Geospatial calculation with GIS background gives full opportunity for new research, especially of the objects in the Solar system. Development of satellite technology for investigation of planets and natural satellites became important after sending the first artificial satellite and human being to the Moon. Other signs of development became visible before the eighties, such as computer methods with the addition of image processing and GIS analysis, which created complete geodatabase. One of the methods of GIS analysis has been used for representing dispersion and a number of phenomena in some areas Bíl et al. (2012); . In this research we used spatial methods for determining dispersion and properties of the craters on the Moon. These methods can be successfully applied along with GIS methods. The methods of sub-pixelization and shadowing may give very interesting and scientific results. With these combined methods we can accelerate the digitalization process.

Special methods for determining craters distribution and properties
After deeper analysis of the craters on the Moon, we created new special methods for determining properties and distribution of the Moon area. Some craters cannot be easily investigated with standard remote sensing techniques. One of the main reasons is the lack of high quality photos and pixelization process. Also, more than 90% of photos, although given in more accurate raster form, are not good enough for absolute digitization. In this paper, we have 1569 craters with different properties and diameter, as well as areas like plates, highlands and mountains. With the finest and precise methods, we have been able to find potentially better places for landing on the Moon. The first method we used is the method of shaded pixelization. We used grayscale blue imagine photos in our investigation. This is the image which presents Lunar Orbiter Laser Altimeter with the resolution of 128 pixels per degree. With the self-created method of shadow, we can get a 1 pixel per 1 second, as well as 3600 times more precise estimation taken from the photos. This special method we called sub-pixelization shade method. Equation for determining percentage into phenomenon of one pixel is where (h ij ) is the distance between the pixel x i and the neighboring pixel x j . We applied a similar soft combined sub-pixel to the topographic maps.
On the (Figure 2). we can see the percentage of some sub-pixels into one pixel.
This method can be called the method of pixel jurisdiction Method of shadow selection can be divided into the main phases and sub-phases. The first phase is transmission and absorption. The number of light rays and the specific wavelength of light can sometimes be absorbed to the high extent, and sometimes partially. For example, if some materials (saturation of pixels) absorb red colors and transmit green and blue wavelength, the combination of green and blue is referred to as a cyan. Reflection occurs when light rays strike a solid object. An object that reflects all light (100%), would be perceived as white. On the other hand, objects which absorb all the light (−100%), are presented black. One of the oldest systems of color properties is the Munsell colour system. This system includes three color dimension: hue, value (lightness), Chroma (color purity). This system was established by Professor Albert H. Munsell in the first decade of the 20th century. Another system is RGB (Red, Green, Blue) system, HSV or HSL, as well as computer graphics opaque ratio. We have one more important color system-CMYK color model. This color model combines three colors-Cyan, Magenta and Yellow. In one bit (1 bit), the number of colors is two, or 2 1 colors. In the standard eight color system (8 bit) or (256 3 ) which produce 16,777,216 different colors. That means that we have more than 16 million colors in standard computer graphic. But in this investigation we have greyscale images of Moon. Those greyscale images have only blue spots which can be used in presenting relief shade. After total analysis we did complete distribution of blue shade on the grayscale background (See Figure 3).
Methods of shadows were used to find better contrast between craters and surfaces on the Moon. The function in this research can be described like function of pixel analysis and its properties referred to as contrast ratio, color saturation, lightness, reflection and distribution of colors. Sub-pixel target mapping was successfully used for soft classified remotely sensed imagery Atkinson (2005). A simple creation of one precise algorithm for sub-pixel target mapping started in the middle of the 80s. All data of land cover must have boundaries Van den Hurk et al. (2003). Even today with modern techniques the accuracy of modern algorithms varies between 80 to 90% Congalton (1991). The accuracy especially increases with new methods of remote sensing like NDVI (Normalized Difference Vegetation Index) or LIDAR (Light Detection and Ranging) data. One of the main limitations is in resolution. Spatial resolution has been the subject of many researches of cover land. Scien-tists have required to evaluate the effect of spatial resolution on delectable spatial variations as characterized by functions of the local variance Woodcock and Strahler (1987). Furthermore, researchers have attempted to find suitable algorithms and methods for better determining of spatial analysis of the pixels Atkinson and Curran (1997), especially in remotely sensed imagery. Early techniques for land cover classification from remotely sensed imagery were focused on hard classification in which each pixel allocated to one class Thomas et al. (1987). Ten years ago and even today, scientists and researchers began to realize soft pixel classification. However, a proportion of pixels will be mixed even where the spatial resolution is fine in relation to the land cover variation (H-resolution case) because some pixels inevitably straddle boundaries between scene objects. The pixel-swapping algorithm is created to receive input of data, and image of land cover proportions k=2 classes. Then, the pixel higher level land proportions are transformed into sub-pixel low and hard cover classes. For example, 20 subpixels by 20 sub-pixels are to be mapped within each pixel, or similarly, the land cover proportion of 77 percent would mean that 77 sub-pixels were allocated to that class. These sub-pixels are allocated randomly within each pixel. Once allocated, only the spatial arrangement of the sub-pixels can vary, not the actual attribute values. Furthermore, the number of sub-pixels allocated within each pixel remains fixed. Firstly, for every sub-pixel the attractiveness of A i of a pixel (i) is predicted as a distance-weighted function of its j=1,2,3,4, 5. . . . . . J, neighbors (see Eq. 1).
Where, Z(x j ) is the binary class of j th pixel at location x j and λ ij is the distance of weight which can be presented (see Eq. 2).
Where, (h ij ) is distance between location (x i ) of pixel (i), for which the attractiveness is desired, the location (x j ), of a neighbor pixel (j), and (a) is a nonlinear parameter of the exponential model. The exponential weightings function is chosen only for color pixels, but this function can be completely used when pixels are in gray scale. The choice of a nonlinear parameter and the number of nearest neighbors are both important. Gray scale pixels are used in the following form (see Eq. 4, Eg. 5-7). The value of gray scale pixel varies between 0 (black) and 255 (white). That could be formulated as 0 1 . . . 255 1 (4) This equation must be presented in short form: Final equation, sustainable for each gray scale pixels is:

GIS analysis and distribution of craters on the Moon
In the earlier stages of the manuscript, we explained how we determined the number and distribution of the craters. By implementing all the craters in the GIS we tried to analyze distribution and properties of craters. After the first analysis, we can conclude where the higher distribution of craters on the Moon surface is (see Figure 4.). In this research, we used QGIS (Quantum Geographical Information System), because this software is robust, open-source and attractive for use. In this software, we also redesigned python code for our own use. Kriging and semi-kriging methods in combination with methods of interpolation may give the best results in geospatial calculations. Although there are a few other methods, priority is given to ordinary kriging and global kriging methods. These methods include autocorrelation or the statistical relationship among the measured points and are very flexible in the presentation of craters' distribution and density. We also used kriging method trough SAGA open source software with the Spatial Analyst Extension Watson (2008); Malczewski (2004); Pew and Larsen (2001).
In addition, in combination with remotely sensed data, a new age of planet cartography can start, especially because new and future spacecraft or satellites will give new results and data. In 20 years' time planet geodesy and cartography will represent new disciplines. We also used geospatial modeling for determining the landing potential of the future spacecraft to the Moon. GIS and modeling of geospatial data are very powerful means for calculating some possible areas, relief characteristics, land anomalies, craters' characteristics, three dimensional profile of the Moon Lew et al. (2016); Tomazos and Butler (2009);Blake et al. (2007); Binns and Nel (2002); Koenig and Abegg (1997). GIS in combination with sub-pixel and pixel analysis can give excellent results. QGIS, SAGA with special tools of geocalculation have been used for representing dispersion and properties of craters. The total number of craters over the whole territory of the Moon was measured by means of remote sensing and sub-pixel method algorithms. By using QGIS within DEM or MDEM (Moon Digital Elevation Model), we also determined zones of altitudes of all the craters.

Color contrast ratio method
For precise determination of craters' position and dispersion we used personally created contrast ratio method and sub-pixel and pixel analysis. Other method which we used is swapping pixel method. These methods were used to determine fine and soft resolution on the maps. The swapping method is very good in the analysis of color, especially for grayscale colors. In the basis of swapping method, we have some equations to explain each method (see .
The coordinates hx and hy present the exponential distancedecay into the window h, T is anisotropy ratio. Anisotropy ratio can be presented in the following form: where ϕ, θ are angels of direction. Transformed coordinates must always be isotropic: where g ′ presents an isotropic model (Thornton et al. 2007). Super-resolution mapping in combination with pixel ad sub-pixel analysis has opened a new horizon in production of finer spatial maps. With the help of pixel swapping analysis we can determine a distance weighted function of its neighbors. See (Eq. 12): Where, W i,j presents attractiveness of sub-pixel (i), to class (j), (M) presents the number of neighboring pixels. Fraction ( j k ) presents the proportion of the neighboring pixel (k) per class (j). Dist. −r i,k presents the inverse distance weighting function between sub-pixel (i). These methods can be applied to determine number, properties and quality of all pixels. In (Eq. 13), we showed mutual relation between pixels: Where (h ij ) presents the distance between the pixel (x i ) and the neighboring pixel (x j ). We combined soft-color contrast method with GIS analysis in our example Valjarević and Živković (2016). Sub-pixel contrast ratio method can be used in color and grayscale analyses. In GIS software, we established the function of historical belt algorithm to find each pixel on the proposed maps. This algorithm presents robust and satisfying method of vectoring raster data and vice-versa . In that way we transformed raster data to vector, in order to find the total number and properties of craters. With the increase of zooming for 12400% and with the contrast ratio of -60 units or with brightness of −160 units we also incerased the maximum of noise on the picture to find border contrast between pixels. In the proposed scale, the total sum of pixels was shown in the area of 1 km 2 . The process of vectorization follows the process of pixelization. After all, in order to find the total number of craters on the Moon, we used several maps of the Moon and Astronomical almanah.

Results and Discussion
The best parts of potential landing on the near side of the Moon were placed in five main belts. These areas are: Montes Jura, Mare Pjato, Terra Pruinare, Mare Frigoris, Tera Fertilitatis. These areas occupy 24% of the total territory of the Moon. Also, we calculated average altitude which would be satisfactory for potential landing. This altitude is between 2,500 and 3,500 m. This belt is 750 m lower than the average altitude on the Moon. Very-low places for potential landing are Mare Homorum, Mare Cognitum, Sinus Iridum, Mare Imbrium, near crater Hiparcus, Tycho, Baco, Ideler, Cambell, Anderson, Mare Crisium, Mare Serentitanis, Mare Tranquillitatis. Low possibility for potential landing are Rimae Hipatia, Sinus Medii, Rima Ariadaeus. The moderate potential belongs to Sinus, Sinus Amoris. Serenitatis, Lacus Bonitatis are places of high possibility. Mare Cognitum and Palus Puterdinis are places of very high potential. The characteristics of very low places show that these places are very close to the edge of the craters and with very high angle of inclination. These angles are between 16 ∘ and 45 ∘ , and are therefore very dangerous for potential landing. The low possibility is within the areas into the craters with diameter between 1,000 and 30,0000 m. Also, angles in these craters have inclination between 10 ∘ and 35 ∘ . Other areas like plateaus with form of calderas, created by lava in this belt, have a great number of fractures. In the moderate areas we have combination of mountains, plateaus, long craters with diameter between 10,000 and 60,000 m. Average inclination of angles is between 8 ∘ and 20 ∘ . Places with very long craters and lava-origin plateaus have high possibility for potential landing. Inclination of angles in these craters is between 5 ∘ and 15 ∘ . Radius of these angles presents good possibility for potential landing, since there are some places, within this belt, with angles lower than 3 ∘ . There are excellent places for potential landing in the north-west and the north-east territories. These territories belong to long-lava plateaus and craters with diameters smaller than 300 m and inclination of angles lower than 3 ∘ . In this belt we have strong possibility for safe landing and future spacecraft expeditions on the Moon. In geological sense the excellent belt for landing is older than others. Reason for that can be the deposition of lava materials in creation of plateaus. Lava filled older craters with materials which can produce alignment of terrain. Absence of massive mountains and deep craters makes this belt extremely valuable and recommendable for future human activity on the Moon. Finally, total results of GIS analysis showed that 12% of the territory is excellent for landing, 18% has very high possibility for landing, 25% moderate, 15% low and 30% very low possibility for potential landing. In the excellent territory average altitude is between 2,000 to 3,000 m. In the territory with very high possibility average altitude is between 3,500 to 4,500 m. In the territory with moderate possibility, average altitude is between 4,500 to 5,2000 m. In the territory with low possibility average altitude is between 5,200 to 6,000 m and in the territory of very low possibility average altitude is higher than 6,000 m. Azimuth in the belt of excellent possibility is 120 ∘ , in the very high one is 230 ∘ , in moderate 150 ∘ , in low belt is 130 ∘ and in very low 270 ∘ . The North-east quadrant and the north-west quadrant are the best with the cover area of 1,5 million of kilometers. This presents only 3,9% of territory of the Moon, but it is enough for landing in the next 1,000 years. The Moon habitat could be built with the help of future precise numerical geodesy and GIS analysis, as well with the new satellite data according to new planned missions.

Conclusion
Micro GIS analysis shows precise geographical and geodesy aspect of the Moon, with the help of online data. Although artificial satellite era and lunar exploration began in the mid-1960s, the Moon and other planets in the Solar system have many long-standing and unanswered questions about planetary and satellite formation. It also refers to the investigation of Lunar origin and geodesy, especially for potential landing. The recent and future development of planetary and lunar geodesy, by means of remote sensing and fine spatial and spectral resolution, will bring about new opportunities to explore and understand the Moon and planets more specifically. Furthermore, with the latest methods and techniques and with free open source data, many scientists can improve their knowledge about the Moon properties. The Moon must be the first station for human spacecraft, therefore, we have set the goal to find the most suitable places for its potential landing. GIS analysis of the Moon, along with astronomical and geodesy data, makes powerful tool to perform this. Private companies and missions will also be included in this project, in an attempt to reach the Moon. With specific cartography data, and new GIS analysis, this job is very likely to have the positive outcome. Joined together, astronomers, geologists, geo-physicians, aeronautic engineers, cartographers, architects and engineers may be included in new expedition after 47 years. The potential new lending with human crew requires the establishment of permanent habitats that have potential for self-expansion and self-sustenance. Moon habitats must contend with surface conditions that include almost no oxygen in the air, extreme cold, low pressure, and high radiation. Human survival on Moon would require living in artificial Moon habitats with complex life-support systems. But before that potential landing must include mission without human crews. Alternatively, the habitat may be placed underground. To contend with these constraints, architects must work to understand the balance between materials and extreme weather on the Moon.

Disclosure Statement: No potential conflict of interest was reported by authors
Data Availability Statement: The [CSV] data used to support the findings of this study are included within the supplementary information file(s).