Digital Navigator on the Seas of the Selden Map of China: Sequential Least-Cost Path Analysis Using Dynamic Wind Data

During the age of sail-powered ships, the maritime trade networks of Southeast Asia were highly cyclical in nature due to the biannually switching wind directions of the East Asian Monsoon. The Selden Map of China provides us with a glimpse of these connections in the early seventeenth century, and it is drawn in a unique way that allows the sailing durations between ports to be measured. In this paper, a novel method of simulating directed sail-powered voyages is developed. The method utilizes ArcGIS Pro’s functionality through Python macros, and unlike the previous least-cost path (LCP) sailing models, it is based on sequential LCP analysis using dynamic real-time series wind data. The optimized routes and sailing durations generated by the macros are then compared against the Selden map. In general, the model performs reasonably well in favourable winds, but is unable to simulate tacking properly in adverse conditions. The results allow the visualization of wind patterns in terms of time spent at sea and demonstrate the inherent natural rhythm of maritime movement and trade in the South China Sea region. The macros are freely available and can be modified to simulate directed sailing in other time periods, localities, and environmental settings.


Introduction
On the 9th of February 2020, I was on the Åland Islands visiting our upcoming teaching excavation site with my colleague. During that day a brisk southerly wind picked up and by the evening it had strengthened into a storm. The ferries travelling from Stockholm via Åland to Helsinki could set out and brave the storm on the open sea but not enter the Mariehamn harbour. We were marooned! In the modern world, it luckily meant arriving in Helsinki only 24 h late, but even that had a ripple effect on the lives of several people because I had to cancel a lecture scheduled for that day.
In the era of sail-powered ships, the same environmental factors were in play and had an even more consequential role in seafaring. In his influential work, Braudel (1995: 355-379) calls distance the first enemy but also points out that the speeds at which those distances could be covered in the premodern world would vary greatly due to weather. One's reach was not constant but fluid, and time more than distance was the traveller's friend or foe depending on the circumstances. This was very much the case in Southeast Asia where the weather is governed by the East Asian monsoon. During the winter monsoon (November to March), winds are northeasterly, but for the summer monsoon (April to September), they blow in the opposite direction. This gave sailing and trade in the area its natural rhythm, where crews might have to wait considerable amounts of time for favourable winds. In addition, short-term changes in local weather had significant effects on sailing. For example, in Yü Yung-ho's account from 1697, thirteen ships set sail at the same time from Xiamen, China, to Tainan, Taiwan, with half of them making the voyage in c. 4 days and the rest arriving in piecemeal during the following ten (Thompson, 1968: 177).
The Selden Map of China from the early seventeenth century offers us a glimpse of the Southeast Asian maritime trade network of that period. The map is drawn in a unique way that makes it possible to calculate the time it took to sail between ports. The trade routes marked on the map show not only the intraregional network, but also how the area was an integral part of the Maritime Silk Route connecting the East and the West. Situated between the economic and cultural powerhouses of China and India, Southeast Asia had a commanding view over the Route by the necessities of geography and wind patterns. Reid (1988: 16-24) has called the period 1570-1630 "the boom years" for the region: it was the producer of highly sought-after luxury products like spices and silk, but also large quantities of bulk goods such as rice, salt, fish, wood and ceramics moved through the networks. Alongside the cargo travelled also people conveying new ideas and ideologies with them.
In this article I propose a novel method of modelling directed sail-powered voyages based on dynamic deterministic wind data in order to nuance the discussion on the early seventeenth century maritime routes and trade in Southeast Asia. It is based on ArcGIS Pro's least-cost path (LCP) functionality but made sequential with the help of Python macros, i.e. small computer programs that automate tasks. The results are compared to the South China Sea routes present on the Selden map ( Figs. 1 and 2). The main aims presented in the article are (1) the development of the method and the macros needed to run the analyses, (2) a summary analysis of outbound routes from Quanzhou to other ports in the South China Sea and a more detailed simulation to both directions on the route Quanzhou-Mekong Island, and (3) the comparison of the simulated results to the sailing durations measured from the Selden map. The focus is solely on the open sea portion of the trade network.
The question of riverside ports situated further inland is beyond the scope of this paper.
In what follows, I focus on studies that share key features with this work: the use of firstly points of sail or polar graphs, and secondly dynamic environmental data. The former are essential in the modelling of sailing speed, and the latter makes it possible to simulate the short-term fluctuations in environmental data that go unnoticed in methods based on longer-term static mean values. This demarcation narrows the field considerably. Points of sail (Arcenas, 2015) and polar graphs (e.g.   Callaghan, 2003a;Davies & Bickler, 2015;Di Piazza et al., 2007;Di Piazza, 2014;Evans, 2008;Fitzpatrick & Callaghan, 2013) depicting the performance of a sailing ship in different winds have been commonly used for the ship speed calculations. The common denominator of these simulations is that they required specially built software.
Although Indruszewski and Barton (2007) used GRASS's wildfire spread routine and ArcView to make anisotropici.e. the cost is dependent on the direction of travel -seafaring models of the Viking Age Baltic sea, Leidwanger (2013Leidwanger ( , 2014see also 2020: 139-151) was the first to introduce the application of ArcGIS's Path Distance tool and its horizontal factor (HF) parameter for this task. His models concentrate on the ancient seafaring in the eastern Mediterranean and are based on modern yearly mean wind direction and speed. The results are presented as time bands, i.e. the sailing time needed to reach every other cell in the research area if an optimal route is chosen. Jarriel (2018) has used a similar approach with monthly mean winds but for paddled vessels without sails. Alberti (2018) has automatized the process further for sail-powered vessels by creating an ArcGIS toolbox. In his case study, Alberti compares Mediterranean historical voyages in classical antiquity to his simulations based on modern monthly mean wind data. In general terms, the modelled results matched reasonably well with the historically recorded durations, with a slight tendency towards overestimation (ibid.). In addition, Safadi and Sturt (2019) have also applied points of sail in the creation of cartograms (see also Di Piazza, 2014) based on sailing durations in the Eastern Mediterranean but their method does not require the use of a HF.
Some models use dynamic environmental data, where wind and/or current speeds and directions change during the voyage affecting the course taken. Levison et al., (1972: 235;1973: 12-17) used a method based on wind roses: each datapoint was divided into 16 compass points and monthly average wind and current speeds and the relative frequencies of their directions were calculated for them. At the start of each simulated day, the values are selected from the table randomly weighted by the frequency distribution. Similar stochastic (i.e. probabilistic) methods have been popular in drift voyage modelling (e.g. Callaghan, 1999Callaghan, , 2001Callaghan, , 2003bCallaghan, , 2011Callaghan, , 2015Callaghan & Bray, 2007;Irwin et al., 1990;Montenegro et al., 2016; see also Davies & Bickler, 2015: 219;Evans, 2008).
In contrast to the stochastic method, a deterministic method based on using realtime series environmental data can also be employed. Montenegro et al., (2006Montenegro et al., ( , 2008 have done this for drift and paddled drift voyages in the Pacific and Atlantic. For the earlier study, the update interval for vessel position and the environmental data is unclear, but the latter uses 2 days for the vessel drift model (Montenegro et al., 2008: 359). In Evans' (2008) simulation of a Polynesian directed voyagei.e. the route end point is known in advance -the wind data is cycled on a six-hourly basis. Perhaps the most interesting parallel to my work in this regard is the modelling tool created by Davies and Bickler (2015). The tool is very versatile allowing many different styles of seafaring, but here I concentrate on its directed sail-powered portion. The following information is based on personal communication with Davies (2021: pers. comm.). The tool uses wind data at a 0.25° resolution, which is updated every 6 h. The ship's position, heading, and velocity are refreshed hourly based on the wind data and are assumed constant for that duration. Winds are typically highly autocorrelated, and consequently the effect of this on the results is small at the same time making the computations considerably faster. The ship speed is calculated from the wind speed and direction using a sailing efficiency function and a polar graph. If a positive speed towards the goal is found, the ship will stay on that heading; otherwise, a tacking direction minimizing the deviation from the direct course is chosen. For sail-powered directed voyages, it has been generally assumed that the navigator can account for leeway and the ocean currents would have a negligible effect on sailing speed, but the tool does allow the inclusion of current data as well if needed. More recently Slayton (2018), in co-operation with Hildenbrand (2015) and Athenstädt (2018), has applied an isochronic method for directed paddle voyages in the Caribbean. In her case study, the environmental data is changed and new isochrones are created every 3 h (Slayton, 2018: 67). In addition, Warnking (2016) and Gal et al. (2021) have applied navigation software designed for modern sailing vessels to model ancient seafaring.
In comparison to the above, the model presented in this paper falls into the categories of directed type of movement and deterministic environmental inputs. It is inspired by the approaches of Leidwanger (2013) and Alberti (2018): the LCPs are calculated using ArcGIS and the points of sail are simulated using the HF functionality of the Path Distance tool. In contrast to the static mean winds used in these studies, dynamic data is employed in a somewhat similar manner to e.g. the Davies and Bickler's (2015) model.

Selden Map of China
The map known today as the Selden Map of China or Dongxiyang Hanghai Tu (Nautical Chart of the Eastern and Western Seas; Fig. 2) is first mentioned in the will of London-based lawyer John Selden in 1653 (Batchelor, 2013: 42;2014: 134;Brook, 2015: 20-22;Nie, 2019: 19-20). After his death the following year, the map has been in the collections of the Bodleian Library (MS.Selden supra 105) since 1659 (ibid.). The map's historical significance to the maritime trade networks of East and South-East Asia was recognized only in 2008 when the trade routes criss-crossing its seas were noticed (Batchelor, 2013). The map remains somewhat of a mystery as its author, purpose, provenance, and exact date are at least partly lost to history. The map was produced in the late sixteenth or early seventeenth century, with scholarly opinions ranging from the Wanli era (1566-1602) to 1624 (Brook, 2015: 173;Nie, 2019: 17 with references). Out of the proposed options, Batchelor's view of a late date c. 1619 based on geopolitical changes present on the map (2013: 55; 2014: 139) seems to be the most plausible. The Selden Map differs from other Chinese maps of similar age in three key ways: (1) in addition to China, it prominently depicts Southeast Asia; (2) there are several sailing routes marked on the map; and (3) it has a scale bar and a compass rose (Batchelor, 2013: 43-49;Brook, 2015: 11-13, 87-88, 105-109;Nie, 2019: 9-14). The map itself is skilfully hand-painted on paper and quite large (c. 158 × 96 cm), so it was more likely meant for use at a prosperous merchant's household (Batchelor, 2013: 42;Brook, 2015: 171-172;Nie, 2019: 9) than as a navigational aid on a ship.
The trade routes present on the map cover mainly an area spanning from Japan in the northeast, to Timor-Leste in southeast and to Aceh, Indonesia, in the west (Batchelor, 2013: 37-40). Additionally, a route connects Aceh to Calicut, India. From thereon, the sailor is guided towards the Arabic world with written instructions to Ormuz, Dhofar and Aden (ibid.). Textual references of a Spanish route crossing the Pacific can be found on the map near Luzon, Philippines (Batchelor, 2013: 45, endnote 17;Brook, 2015: 121-122). Quanzhou, China, has been identified as a central link among the ports portrayed on the map (Batchelor, 2013: 49;2014: 51-52;Brook, 2015: 113;Nie, 2019: 32-33). It is also noteworthy that not all ports marked on the map are connected to the network suggesting a distinction between short-sea and long-distance shipping (Batchelor, 2013: 49). Many of the routes do follow coastlines, but some also brave open seas. The routes are made up of straight segments, most of which are accompanied by a label declaring a compass heading. These labels correspond to the notation used in the map's compass rose (Batchelor, 2013: 47-49;Nie, 2019: 10-11).
The compass rose and the scale bar are situated at the top of the map. The compass rose has two sets of compass points and their labels marked on it: one is divided into cardinal and intercardinal directions, and the other into 24 compass directions in a fashion typical to Chinese compasses (Batchelor, 2013: 46-47;Brook, 2015: 87-88;Davies, S., 2013: 97-98;Nie, 2019: 14). The scale bar consists of ten major divisions, which are further partitioned into ten subdivisions (Batchelor, 2013: 46). The unit of measure has not been marked on the map, but according to the current consensus first proposed by Davies (2013: 98-99), the minor divisions depict 2.4-h ship's watches (geng) and the major refer hence to days of travel (Batchelor, 2013: 49;Brook, 2015: 162-163;Nie, 2019: 14-15). Davies (2013: 99) bases his viewpoint on Chinese route guides using the same compass direction/geng system and the lengths of the route segments falling within the typical performance of traditional junks. Batchelor (2013: 49) also points out that the written instructions from Calicut indicate the use of geng on the map. Hsu (1988: 102) has described geng aptly: Gen[g] was a time unit in historic China; one day and night constitute ten gen[g], or one gen[g] equals 2,4 hours. In navigation, the expression 'X gen[g]' or watches was used to mean that the X time was needed to travel the Y miles between place A and B. Therefore gen[g] is an expression of time/ distance. Since ship velocity was affected by several factors, it is impossible to define the gen[g]-mile equivalence precisely.
If understood as a distance unit, the estimates range between 14 and 20 miles or 22.5-32 km (Hsu, 1988: endnote 2). According to Davies's calculations, the speeds achieved in different parts of the Selden map vary from one to nine knots in conjunction with the typical winds for the area (Davies, 2013: 99). He provides us a single example from the mouth of Manila Bay, Philippines, to Brunei: the distance of c. 1185 km and the duration of 6.2 days represents an average speed of 4.3 knots. In comparison, the average speed for all routes calculated from the rutter Shunfeng Xiangsong (Fair Winds for Escort, also known as the Laud rutter; MS.Laud Or.145) is 4 knots (Mills, 1979: 72). It dates to the early fifteenth century and covers partly the same routes as the Selden map (ibid.).
Both the compass rose and the scale bar are tilted c. 6° counterclockwise in comparison to the edges of the map. Contrastingly, an empty bordered rectangle with the same orientation as the paper has been drawn next to them. These features suggest it to be the mapmaker's way of designating magnetic declination (Davies, 2013; see also Batchelor, 2016), and the route lines seem to follow the same tilt (Brook, 2015: 164-165). In general, most of the lines coincide with their stated compass headings within 7.5° margin of error (Kogou et al., 2016: 14). The mapmaker achieved this accuracy by drawing most of the routes first and adding the landmasses later around them (Batchelor, 2013: 37, 41;Brook, 2015: 161-162). This became evident when a draft of the route stretching from the Goto Islands, Japan, along the Chinese and Vietnamese coasts to Pahang, Malaysia, on the otherwise almost empty verso of the map was found during conservation (Minte et al., 2014: 117), and subsequent analysis done on the pigments on the recto has confirmed this to be the case for the majority of the routes (Kogou et al., 2016: 12-19). The largest deviation can be found on the aforementioned Aceh-Calicut route, most likely due to the mapmaker running out of space at the edge of the paper (Kogou et al., 2016: 14). Otherwise, the map is rather accurate for its time (Brook, 2015: 166, image 26).
The Selden Map of China represents in many respects major advancements in Chinese mapmaking, and it is likely that some of its features such as the compass rose and scale bar were influences from Western cartography (Batchelor, 2013: 43;Brook, 2015: 109;Nie, 2019: 12-14; see also Kogou et al., 2016: 20 for possible Western traits in the materials used to make the map). The mapmaker considered the routes to be an integral feature of the map (Batchelor, 2013: 47-49;Davies, 2013: 100). They are not just rough graphical illustrations; rather, they function as a visual rutter and represent actual sailing routes (ibid.). Presumably, the route lengths depict typical sailing durations in favourable winds. The map was most likely a commercial navigation chart owned by a merchant or a trading organization (Batchelor, 2013: 42, 56-57;Brook, 2015: 167) and should therefore probably be understood as a representation of the mercantile interests of a specific party rather than a complete description of all trading routes of the time. In the vast majority of cases, the route lines do not extend to the ports or the shoreline but begin some distance offshore. According to my interpretation, the sailing instructions would guide the mariners to suitable areas nearby ports, from where they would use other means of navigation to reach the ports themselves. This rationale is exercised throughout this article: when distances and routes are calculated, their start and end points are slightly seaward as they are on the Selden Map.
This short description of the Selden map is based on sources in English and covers only the aspects relevant to this article; for example, discussions on the authorship and provenance have been omitted. For a more comprehensive overview, see the works of Batchelor (2013Batchelor ( , 2014Batchelor ( , 2015Batchelor ( , 2016, Brook (2015), and Nie (2019). In the following sections, I construct a sequential LCP sailing model for the South China Sea and compare its simulated voyage durations to the values measured from the Selden map. The aim is to gain a clearer understanding of seafaring in terms of time, and its implications to the interconnectivities in the region.

Method
The modelling methods mentioned in "Research History" have all their own pros and cons. In short, stochastic environmental data can account for the variability found within, but random selection does not follow the strong temporal autocorrelation apparent in meteorological processes (Levison et al., 1973: 16). On the other hand, long-term annual or monthly means give us a generalized picture of the actual weather patterns but average out shorter fluctuations, storms, and other random events in the data. These again are present in deterministic data, but the simulations using this type of data usually cover relatively short time periods, which might not reflect the typical circumstances. Approaches similar to the directed sailing portion of Davies and Bickler's (2015) tool are good at simulating known courses, but do not try to find the optimal route. Methods based on finding LCPs accomplish this, but the optimized paths might not reflect the actual courses taken by the mariners of the past.
The method used in this paper is based on cost distance and LCP analysis functionality of ArcGIS Pro. In cost distance analysis, the research area is divided into a regular grid and a cost (or a NoData value if a cell is to be excluded from the analysis) is assigned to each of its cells. The cost can be for example energy consumption, price or -as in our case -time spent moving across a cell. In its basic form, cost distance calculates for each cell the minimum accumulated cost needed to reach it from the source point (Fig. 3). In ArcGIS's LCPs, movement from one cell can happen only to one of its eight neighbouring cellsi.e. to cardinal or intercardinal directions -and the longer distance covered in diagonal movement is taken into account in the calculations (Esri, 2019a. Though computationally easy, this approach does suffer from grid-induced biases: the accumulated cost bands radiating out from the source point are not perfectly shaped but tend towards a somewhat polygonal shape (e.g. Dunn, 2010: fig. 1;Etherington, 2016: 48-49), and straight lines running in other angles than perpendicular or diagonal to the grid cannot be perfectly portrayed (e.g. Antikainen, 2013: fig. 1; Herzog, 2013: 191-192).
ArcGIS Pro's Path Distance tool adds to this the option of using vertical and horizontal factors to make the analysis anisotropic (Esri, 2019b(Esri, , 2020c(Esri, , 2020d. The vertical factor is typically used in land-based analyses to simulate the effects of slope on movement: in general terms, it takes more effort or time to travel uphill than downhill (for recent studies see e.g. Alberti et al., 2020;Hart et al., 2019;Seifried & Gardner, 2019). On the other hand, the HF can be used to introduce horizontal frictional elements: e.g. a car will experience more drag travelling upwind than downwind. In our case, no vertical factor is needed as the sea can be assumed to be flat, but the HF in combination with wind direction data is used to scale the ship speed according to the points of sail (Alberti, 2018: 512, 515;Leidwanger, 2013: 3305;Fig. 4). Because only HF needs to be considered, the cost distance between neighbouring cells a and b is calculated in the following manner (c.f. Esri, 2020d), where CD is cost distance, CS cost surface, HF horizontal factor, and SD surface distance: The HF can be defined in the Path Distance tool as a table for the angles between 0 and 180° ( Fig. 7), the values of which are then mirrored to cover the other half of the circle (Esri, 2019b(Esri, , 2020c. A factor of one leaves the cost unchanged, while higher values increase it and lower values have the opposite effect (ibid.). The Path Distance tool produces two new grids: a cost distance raster as discussed above and a backlink raster (Esri, 2020d. The backlink raster's cells record the direction to a neighbouring cell along the path with the least accumulative cost towards the source point (ibid.). The Cost Path tool then extracts this LCP between the source and destination cells (Esri, 2020a(Esri, , 2020b. During the Path Distance process, the wind speed and direction rasters cannot be changed. In order to include dynamic winds to the analysis, I have programmed a macro that performs it in several 6-hourly segments. This duration was chosen, because it is the measurement interval in the wind dataset (see "Environmental Data") and in order to keep the computing time feasible. In short, the macro functions as follows: (1) path distance and LCP are calculated; (2) the segment with accumulated cost less than or equal to 6 h is selected from the LCP; (3) the LCP cell in this segment with the highest accumulated costi.e. the furthest point reached on the LCP within 6 h -is chosen to be the new source cell; (4) the wind data is   Fig. 3 The principles of cost distance and LCP analyses (Esri, 2019a(Esri, , 2020d. The effects of a horizontal factor are omitted and only the isotropic cost distances are shown here for the sake of clarity advanced by 6 h; and (5) the process is repeated until the destination cell is reached (for further details, see "Macro 2: Digital Navigator"). In addition, a new ship can be launched from the source cell at a selected time interval, for example daily.
The major drawback of this method is the difficulty of incorporating ocean currents into the analysis. Leidwanger (2013: 3304) has argued that if currents are generally slow, they only have a minor impact on ship speed compared to the wind. Alberti (2018: 513 with references; see also Di Piazza et al., 2007: 1221Newhard et al., 2014: 383) shares the same view and adds that superficial currents are in general wind-induced and tend to move into the same direction as the wind. The South China Sea is a semi-enclosed sea, where currents are mainly wind driven (Potemra & Qu, 2010: 117), and the winter and summer mean surface currents are mostly slower than 0.6 knots on average (Park & Choi, 2017: fig. 2). On the grand scale, the upper ocean circulates in the Southeast Asian Seas as follows: the water flows from the Pacific to the South China, Celebes, and Halmahera Seas; makes its way from there to the Java, Flores, and Banda Seas; and exits into the Indian Ocean mostly via the Lombok and Ombai Straits and the Timor Passage (Potemra & Qu, 2010: 118-120). Overall, the main inflow sources are the Mindanao and New Guinea Coastal Currents (ibid.). In addition, the South China Sea also receives water from the Kurishio Current through the Luzon Strait (ibid.). The Luzon Strait inflow is at its height during the winter months creating -with the help of the northeasterly monsoon winds -a current flowing down the coasts of China and Vietnam and large-scale counter-clockwise circulation of the South China Sea (ibid.). On average, the speed of the coastal current is c. 0.3-0.6 knots (Park & Choi, 2017: fig. 2). Furthermore, the wave direction generally follows the wind direction and wave heights are less than 2.5 m for the majority of time (Osinowo et al., 2017). Based on these traits, it was deemed acceptable to exclude currents from this model.
Turning now onto the sailing properties of the simulated ship, it should be emphasized that properties used in this study are simplifications and they are not meant to represent the characteristics of a specific wreck or vessel; rather, they function as rough generalization of the performance of early modern period sailing ships in Southeast Asia region to form a basis for discussion. For this simulation, the ship's velocity is mainly derived from wind speed and direction by its ship-speed-to-windspeed ratio and the sailing performance in different points of sail. Other factors such as hull form are not taken into account.
In general terms, a sail-powered ship is propelled forward by the wind pushing against its sails or by the sails acting like a wing and generating lift. It naturally follows that good speeds can be achieved sailing downwind. Conversely, only modest gains upwind can be made by beating to windward, i.e. tacking between closehauled points of sail in a zig-zag course. Making headway directly against the wind -into the aptly named no-go zone -is impossible. Counterintuitively, a heading parallel to the true wind direction is not the most efficient point of sail: on a running course, the boat is also running away from the wind and the effective apparent wind felt onboard is lower. Instead, the boat will move the fastest in broad reach with the wind blowing at an angle behind the ship or in beam reach when the wind is perpendicular to the course. The angles and speeds achieved in the different points of sail are dependent on the hull shape and sails of the ship and therefore vary from vessel to vessel.
Many of the ships carrying cargo in Southeast Asia in the early seventeenth century were Chinese junks or Southeast Asian jongs (e.g. Manguin, 1980;Manguin, 1993;Reid, 1988: 18-19, 36-43). Although polar graphs are available through experimental archaeology for other ship and sail types (e.g. Brandt & Hochkirch, 1995: fig. 6-9;Nomoto et al., 2003: fig. 13), I was unable to find numerical performance values for historical junk rigs during my research and therefore a comparable substitute had to be selected. First, let us consider shortly what can be said on the capabilities of the junk rig. For example, Du Halde (1738: 329), basing his view on Jesuit missionaries' reports, assessed the junks thusly: [...] they proceed very swiftly, when steering right before the Wind ; nay, can, if we believe them, keep up with our best Sailors, and even leave them behind. But then with a quarterly or Side-Wind they cannot hold it, and are driven out of their Course : not to mention the Danger they are in of being turn'd about, when they are surpriz'd with a sudden Flurry of Wind.
His description is mostly confirmed by a reconstruction of an Early Qing Dynasty (1644-1799) Ganzeng junk of war, the Princess Taiping, which managed to reach average speed of 6 knots (Lu & Hao, 2011: 8). The ship was reported to be stable but had a significant transverse drift typical to most Chinese junks (ibid.). In sea trials with a one masted modern dynamometer boat equipped with a Chinese lugsail, the tacking manoeuver was very easy and could be done c. 52.5° from the wind (Masuyama et al., 2005: 128).
Whitewright has analysed historical accounts of voyages made under square and lateen/settee sails and offers potential point of sail performance for both (2011: 13-15, fig. 5; see also Casson, 1971: 281-296). Whitewright utilizes the concept of velocity made good (VMG), which he defines as the absolute speed of a vessel over a direct course between two points (2011: 3-4). For example, if a ship is tacking directly windward, its speed towards the tacking angle can be higher but the VMG denotes the progress made towards the destination. The results for the square and lateen/settee sails are almost identical: on reaching and running courses, the VMG is 6 knots on both; in broad reach, the maximum VMG for square sail is 12 + vs. the lateen/settee's 10 + knots and the square sail is slightly faster at a close-hauled course (< 2 vs. < 1.9 knots) but lateen/settee can at its best sail slightly closer to the wind (c. 56-67° vs. c. 60-65°) (ibid. : 13-14, fig. 5). The average VMG for both sail types has also remained essentially the same throughout the last two millennia (ibid.: 14-15, fig. 6); as Englert (2006: 177) has pointed out, on a running course, a 900-year-old Viking design can still sail neck and neck with a modern Bermuda-rig racer without its spinnaker! Out of these two options, the lateen/settee sail with 67° upper limit on close-hauled was selected for the simulated ship as the more conservative options (Fig. 4). Although Whitewright focuses mostly on the Mediterranean, his dataset includes also Arabic fifteenth century lateen/settee sail voyages in the Red Sea (Whitewright, 2011: 11-13).
To compare the effects of lateen/settee sail close-hauled upper limit of 67° versus of 56° from the wind direction, the sailing durations based on the winds of April 1979 (see "Environmental Data" for the reasoning behind choosing this year for analysis) were simulated for both values with otherwise identical settings on the route Quanzhou-Mekong Island ( Fig. 5; see "Macro" for simulation details). In favourable conditions -which are the main concern of this article -from the 1st to the 20th of April, the difference is on average 0.2 (0.2) days (see "Macro 2: Digital Navigator" on why two values are presented). In some cases, the less capable ship is even slightly faster most likely due to the superior ship's digital navigator venturing more eagerly to adverse winds. On the other hand, the difference becomes apparent with the unfavourable winds later in the month when the gap grows to 15.2 (16.6) days on average. On more complicated routes, the effect is likely to be somewhat more pronounced.  (1995), Harries et al. (2000), and Nomoto et al. (2003) that the shipspeed-to-wind-speed ratios of traditional vessels fall typically to the range of 0.25 to 0.5. For example, the ratio for Naniwa-maru, a replica of an eighteenth to midnineteenth century Japanese Bezai type trader, was c. 0.3 in a real-life testing run. At the time, the ship's bottom was fully fouled with 1 cm barnacle layer and the authors conclude that the ratio would have been 0.4 in broad reach were the hull clean (Nomoto et al., 2003: 143-144, fig. 13). Palmer (2009a: 65) has also pointed out that the ratios of vessels slow for their size -such as traders -are likely to be reasonably constant and tests done in one set of circumstances can be extrapolated to different wind speeds. Therefore, a mid-range constant ratio of 0.4 was chosen for this simulation.
Previous LCP sailing models based on ArcGIS's Path Distance tool have used long-term yearly (Leidwanger, 2013: 3304) or monthly (Alberti, 2018: 519) mean winds where storms and lulls are averaged out and hence no extreme winds will be encountered. In this study, actual measured winds from specific moments in time are used and our digital navigator may find himself in a tempest. Sailing becomes more strenuous in high winds, and therefore a method lowering the ship's VMG in relation to the increasing wind speed in these circumstances must be defined. For example Arnaud (2005: 19-20), writing from an ancient Mediterranean perspective, states that making headway even downwind becomes difficult in winds greater than force 5 on the Beaufort scale, and from force 7 onwards the situation quickly becomes critical. Alberti's (2018: 514, fig. 3) rescaling method based on the maximum wind speed cannot be applied here -the 6-hourly changing winds and the possibility of storms would cause the ship speed at a specific wind speed to fluctuate -and therefore another way of tackling this problem needs to be devised. Armed with the information above, the ship's VMG in broad reach is set to increase linearly until peaking at 25 knot wind (mid force 6), after which it starts to suffer from the gale and heavy seas declining again linearly up to 40 knot wind (high force 8). From thereon, the VMG is set to a constant low value to allow modest gains and ensure the ship's arrival to destination at some point even in adverse conditions. For the other points of sail, the values are calculated from the broad reach VMG using a HF (Fig. 7). The ship's point of sail VMG in relation to wind speed is presented in Fig. 6.
In summary, given the LCP methodology and the sailing properties used in this model, our digital navigator can be thought to possess certain characteristics: • Has limited all-knowingness. Always knows its position and is free to choose the optimal route in the currently prevailing winds but has to update its plan in step with the 6-hourly wind changes. This allows the navigator to be sometimes caught in unfavourable winds therefore introducing a slight element of randomness and indecision to the route • Sails round the clock at the maximum performance of the ship. Does not need rest, repair or supplies. The skill level of the crew, unoptimal sailing etc. are somewhat represented due the point of sail performance being based on actual historical voyages • Sets sail in all weather conditions but aims to evade extreme winds on the route • Tries to avoid shallows and certain island groups (see "Macro 2: Digital Navigator")

Environmental Data
As historical climatological datasets with high temporal and spatial resolution do not exist, using modern data is a typical practice when modelling ancient seafaring. The digital environmental datasets used in this study come from two sources: Climate Forecast System Reanalysis (CFSR) and General Bathymetric Chart of the Oceans (GEBCO). The CFSR dataset (Saha et al., 2010) is available from January 1979 to March 2011 and its successor CFSv2 (Saha et al., 2014) from April 2011 to near present with hourly temporal and a global spatial resolution of 0.5° × 0.5° (c. 56 km at the equator). The high temporal resolution is achieved by forecast models, which are reinitialized every 6 h (0000, 0600, 1200, 1800 UTC). From the numerous variables available in the CFSR, 10 m above ground u-and v-components of wind (m/s) are used in the simulation to calculate wind speeds and directions. The data is based on satellite measurements (Saha et al., 2010(Saha et al., : 1021(Saha et al., -1022 and represents the near-surface wind conditions. In turn, the General Bathymetric Chart of the Oceans' GEBCO_2020 is a global gridded ocean and land terrain model (GEBCO Compilation Group, 2020). Its spatial resolution is 15 arc seconds (c. 460 m at the equator) with heights/depths stated in 1-m increments from the modern mean sea level. In the simulation, rasters derived from GEBCO_2020 are utilized to differentiate land from sea and to introduce extra movement cost in proximity to island groups and in shallow parts of the ocean. The processing steps of these datasets are described in detail in "Macro". The question though always remains the same: how well does the data depict the circumstances of the past? Modern weather station data from mainland China shows that winds have been weakening in the latter part of the twentieth century: annual mean wind speeds have declined by c. 30% and the number of windy days (mean wind speed > 5 m/s) has dropped c. 60% from 42 days in 1969 to 17 days in 2000 (Xu et al., 2006). Further studies (Fu et al., 2011;Guo et al., 2011;Jiang et al., 2010) have since confirmed and expanded upon this trend. The effect is significant on both East Asian summer and winter monsoons (Xu et al., 2006). Annual wind speed change correlates with latitude: the decline is smaller in southern China (− 0.08 m/s per decade) and greater in the northern parts of the country (− 0.29 m/s per decade) (ibid.). Warming over the ocean and cooling in south-central China contribute to the weakening of the summer monsoon while the winter monsoon is affected by global warming (ibid.). The reasons behind these changes appear to be anthropogenetic, air pollution, and greenhouse gases respectively (ibid.).
Winter monsoon wind velocity (WMV) can be reconstructed further back into the nineteenth and early twentieth centuries by using Sr/Ca ratio (Liu et al., 2008) or oxygen isotope δ18O (Song et al., 2012) (Song et al., 2012). This however is the temporal limit of our current knowledge of wind speeds in the area as monsoonal wind intensities are difficult to infer from geological records (Sun et al., 2015: 133) and wind speed data is not included in Chinese historical records (Wang & Zhang, 1988). In China 1969-2000, the mean wind speeds had a negative correlation with annual mean temperature (y = − 0.335x + 5.515 [R 2 = 0.49, p < 0.001]) and a positive with ground station measured annual mean solar radiation (y = 0.322x − 2.139 [R 2 = 0.51, p < 0.001]) (Xu et al., 2006). The global temperature anomaly of the early decades of seventeenth century compared to the base period 1951-1980 is c. − 0.5 °C (2° Institute, n.d. a; based on Moberg et al., 2005). Regional data is also available for Southeast China, c. 0.0 to − 0.1 °C, but with a base period of 1901-1950 and values missing for late twentieth century (Ge et al., 2010: fig. 3) and thus not directly comparable with the aforementioned values. Moreover, historical values for mean solar radiation reaching the ground are not known, but the total solar irradiance (TSI) hitting the Earth's higher atmosphere in 1610-1620 fluctuated between c. 1360.5-1361.5 W/m 2 (Kopp et al., 2016a: 2958-2960). During the last 50 years, the global mean temperature has risen sharply by a degree and the wind data should therefore be chosen from an early point of time as possible. In the case of CSFR, this is 1979-1980, during which the temperature anomaly was c. + 0.2 to + 0.3 °C (2° Institute, n.d. a) and the TSI c. 1361.6 W/m 2 (Kopp et al., 2016a: 2958-2960). If the correlation between mean wind speed and temperature holds true for the early seventeenth century as well, the c. 0.8 °C difference in temperature would represent c. 0.3 m/s higher mean wind speeds then than in 1979-1980 with the TSI possibly compensating slightly. It is also important to note that 1979-1980 are neutral El Ninõ-Southern Oscillation (ENSO) years (Song et al., 2012: 136, table 2). ENSO is a recurring, though somewhat irregular, climate phenomenon in the tropical Pacific. The periodic warming and cooling of the sea surface temperatures in the central and eastern parts of the ocean cause the so-called El Niños and La Niñas, which have major effects on global atmospheric circulation and therefore affect the wind patterns in Southeast Asia as well.
On the other hand, the mean sea level has been more stable: in the first half of seventeenth century, it was c. 20 cm and in 1979-1980 c. 10 cm lower than in 2020 (2° Institute,n.d. b; based on Kopp, R. E. et al., 2016). Because the z resolution of the GEBCO_2020 is 1 m, the modern shoreline and depth information was used as is. All things considered, the modern datasets chosen for this study have a high enough spatial and temporal resolution for modelling purposes and they are reasonably well comparable to the conditions of the Selden map timeframe.

Macro
The two macros used to prepare the wind data (macro 1) and simulate sailing durations (macro 2) were written in Python for ArcGIS Pro using Spyder development environment (see Fig. 8 for flowcharts). They use ArcPy site package and a Spatial Analyst licence to access the geoprocessing tools of ArcGIS Pro and Arrow library for timecoding. Both macros were run on ArcGIS Pro 2.5 (Python 3.6.9) as 2.7's Path Distance and Distance Accumulation tools had a program crashing bug when using HF tables at the time of writing. At the moment ArcGIS and its users are transitioning to ArcGIS Pro, but the older ArcGIS Desktop is still widely in use. To make the macro 2 also adaptable to the latter, I chose the legacy Path Distance and Cost Path tools found in both over the ArcGIS Pro's new Distance Accumulation and Optimal Path as Raster tools. "Macro 1: Wind Data Batch Processor" and "Macro 2: Digital Navigator" describe how the macros were used to calculate the optimal sailing durations and routes for ships launched each day of the year between Quanzhou and Mekong Island. The macros and a test dataset are available for download at Zenodo (Perttola, 2021a(Perttola, , 2021b.

Macro 1: Wind Data Batch Processor
The macro 1 processes the CFSR wind data into wind speed and direction rasters. Before running, the hourly CFSR u-and v-components of 10 m above ground wind (m/s) were downloaded as NetCDF files (Asia-Pacific Data-Research Center, n.d.), and the GEBCO_2020 ocean and land terrain model as GeoTIFF (GEBCO Compilation Group, 2020). The GEBCO_2020 was further processed beforehand by resampling it into 500-m cell size while changing the coordinate system from WGS84 to Asia South Equidistant Conic (ESRI:102,029), and reclassifying it into a land/sea raster (land: z ≥ 0 → 1, sea: z < 0 → 0).
When the macro is run, a start time is set and the following actions are performed in a six-hourly for loop. First the u-and v-components are processed separately: vector point layers are created from the NetCDF files, the layers are projected into  Fig. 8 Simplified flowcharts of the macros Asia South Equidistant Conic, the points are interpolated to a raster by using natural neighbour interpolation (Sibson, 1981) with the land/sea raster as a snap raster, and the rasters are clipped into the study area extent. In the next phase, the u-and v-component rasters are used to calculate wind direction and speed. The former is accomplished by the Raster Calculator formula (180∕π) * ATan (u, v) . The outcome is on a scale of − 180-180° and it is converted to 0-360° by adding 360 to the negative values. The latter is calculated with sqrt(u 2 + v 2 ) . Finally, the land areas are set to null on both rasters using the land/sea raster. They are saved to disk as GeoTIFFs and the time is advanced by 6 h before returning to the start of the loop. The loop can be set to repeat the wanted number of times.
The end results are six-hourly (0000, 0600, 1200, 1800 UTC) wind direction (0-360°) and speed (m/s) rasters with 500-m cell size. In this study, the time period 1. 1.1979-1.7.1980 was processed in this manner. In the following sections, these files are used to simulate voyages departing daily between 1.1. -31.12.1979. The wind data for the first half of 1980 is needed for the ships setting sail near the end of the year to complete their voyages.

Macro 2: Digital Navigator
The macro 2 simulates the sailing duration from point A to point B by finding the optimal route using LCP analysis. During the processing, the wind data is changed every 6 h, and the analysis is repeated for ships launching each day of the year. The macro begins by setting the ship characteristics (see "Method" for the values used in this study), the route start and end points, and the start time. Next, a for statement is used to loop the days of the year and reset the route start point to its original value.
The six-hourly sailing duration and route calculation is done within the following while statement that loops until the ship reaches the route end point. In the while loop, the wind speed for that specific hour is rescaled into ship speed in broad reach (Fig. 6). The ship speed is then used to calculate the cost, which is defined as the time used to cross a map unit, in our case seconds per meter. At this point, two extra cost factors are introduced to mimic the navigator's avoidance of these areas: shallow parts of the sea and the proximity to the major island groups mentioned on the Selden map (Fig. 11). Both are derived from the resampled GEBCO_2020 dataset and prepared beforehand by the user. Shallows are defined as the depth of − 3 ≥ z < 0 m and carry an additional cost factor of 1.5. On the Selden map, the Paracel and Spratly Islands are marked east off the coast of modern-day Vietnam (Batchelor, 2013: 52, fig. 9) with the accompanying texts "shoals for ten thousand li in the shape of a ship's sail" and "reefs for ten thousand li" (Brook, 2015: 16-17).
Li is a Chinese distance unit of 500 m and the expression "ten thousand li" carries the meaning "a lot" (ibid.). Brook has interpreted the texts as a warning for the navigator to steer clear of the area (ibid.), and following this view, an extra cost was attributed for the movement in the vicinity of the island groups. The island cost factor is formed by first calculating Euclidian distance up to 50 km from the Paracel and Spratly Islands and using formula 2 − distance∕50000 . That is, the factor is 2.0 when the distance to an island is zero and it decreases linearly down to 1.0 when the distance is 50 km or greater. The total cost in broad reach is reached by multiplying the wind speed-derived cost with the shallows and island proximity factors. In other words, these areas can still be sailed through but at slower speeds and higher time costs.
The cost is then used as an input for the Path Distance tool alongside the wind direction, HF table, and the current start point. It is at this stage that the different points of sail (Fig. 4) and their effect on the cost in the prevailing wind direction are taken into account with the HF (Fig. 7). In this study, the HFs for the different points of sail are calculated using formula HF = (pointofsailmaxspeed∕referencespeed) −1 . The choice of the reference speed is dependent on how the wind speed data is rescaled into ship speed, and as mentioned above, this is done in the macro (Fig. 8) by rescaling the wind speed into ship speed in broad reach. Therefore, the reference speed in this case is the maximum speed in broad reach, resulting in a HF of 1 for broad reach and the HFs for the other points of sail being calculated accordingly. If the rescaling is done differently, other reference speeds are possible (cf. Alberti, 2018: 515-517 where the speed for running and beam reach is used). The no-go HF was subjectively chosen a value of 20 to reflect a very high cost of travel when in irons. Two rasters are created by the tool: a path distance raster, which identifies the least accumulative cost distance for each cell, and a backlink raster that defines the direction from each cell to its surrounding neighbour along the least accumulative cost path leading to the start point (Fig. 3).
The path distance, backlink, and route end point rasters are passed on to Cost Path tool. It calculates the optimal route with the least accumulative cost between the current start point and the end point. This path is then overlaid with the path distance raster and the segment with accumulated cost less than or equal to 6 h is selected. The maximum value of this segment is added to the sailing duration so far, and the cell is selected to be the new start point for the next loop. A small time skip happens at this point: the accumulated cost is usually somewhat less than 6 h, but in the next loop, the winds are advanced by exactly 6 h out of necessity. Therefore, the sailing duration for the full route is stored and presented in this article both as the sum of the segment accumulated costs and the full 6-h loops required to complete the voyage (Figs. 5,9,and 12) with the true value likely being somewhere in between. The segment is further overlaid by wind speed, and its minimum and maximum values on the segment are calculated and they are stored if they exceed the previous values. Finally, the segment is combined to a mosaic raster forming the full route piece by piece.
At the end of the while loop, time is advanced by 6 h and the segment end point is compared to the route end point. If the ship has not arrived to the end point, another loop begins. Otherwise, the end results of the macro -the durations and the minimum and maximum winds encountered on the route and the route itself -are saved into a text file and a raster respectively. After this, the start time is advanced to the next launch day and the for loop launches a new ship from the original start point to sail once more. 1

Results and Discussion
As discussed previously, the Selden map is constructed in a way that allows sailing durations to be calculated from its route lengths. The calculations can be easily and accurately made with the help of GIS. After scaling the Selden map according to its scale bar and tracing the routes as polylines, it is possible to sum up the durations for any route combinations. In this study, the measurements were made from the Bodleian Libraries (MS. Selden Supra 105, Map recto) photo as straight segments, the length of which was rounded to the nearest half geng (cf. Mills, 1979: 84) before summation. The results can vary by a few geng as the exact segment start and end points are slightly ambiguous and there are damaged parts of the map that require some interpretation from the measurer. Additionally, if the route is replicated on a modern map, its geodesic length can be used to determine the average speed on the Selden map. For example, the route from Quanzhou to Mekong Island is illustrated on the Selden map to take 10.85 days of sailing and a similar route on a modern map is c. 2300 km in length resulting in an average speed of 4.8 knots. The modern lengths and consequently the average speeds need to be considered approximations as plotting the Selden map routes on to a modern map is a somewhat subjective task. Also, as discussed in "Research History", the start and end points are marked as is on the Selden map, i.e. some distance away from the ports out to sea. If one were to consider portto-port travel, some additional time would need to be added to the values presented in this paper for the travel to/from the port and docking/undocking procedures.
The Selden map does not provide us with any clues concerning the time of year its each sailing duration is based on. To gain a general view, sailing durations for ships launched on the first day of each month were therefore simulated to several ports in the research area based on the winds of 1979-1980 (Table 1). In six cases out of eleven, the simulated optimized routes underestimate the Selden map durations by a few days, and it is probable that somewhat lower values would be found were the full year or the fastest months processed daily. The Selden map routes are unlikely to be based on the fastest journeys on record at the time but rather on typical speeds achieved in good conditions, and thus some underestimation is to be expected. Other possible reasons for the underestimation are as follows. Firstly, the digital navigator is free to choose the most optimal courses, some of which cut through the open sea unlike on the Selden map, though the coastal route southwest is only marginally slower at least on the fastest month of November (Fig. 5). Secondly, years are 1 A moderately powerful PC (AMD Ryzen 7 3700X, Samsung 1 TB 970 EVO Plus SSD, 32 GB DDR4 3600 MHz) took c. 55 s to process one 6-h cycle on the route Quanzhou-Mekong Island, i.e. 10 days of sailing required c. 37 min to process. This is due to the considerable file sizes and the Arcpy tools used not being efficient in parallel processing. Also, as a recent newcomer to Python and Arcpy, the macros could most likely be further optimized. Altogether, the processing for this article was split between six computers. The processing could possibly be made faster by modifying the macro to run multiple days in parallel on a supercomputer, but this is yet to be attempted. not brothers, and the winds of 1979-1980 could be more favourable than normal. Thirdly, the ship properties used in the simulation could be too optimistic. The largest outlier is the route Quanzhou-Kedah, where the fastest simulated duration of 32.6 (34.0) days overestimates the Selden map value by c. 10 days. Although taking a different route through the South China Sea, the ship passes Johor on day 16, i.e. in the same time as its Quanzhou-Johor sister launched on the same day completes its voyage. The reason for the overestimation is the unfavourable winds in the Strait of Malacca. In general terms, it is apparent that on the analysed year the launch dates in January and October-December would have been the most beneficial for southwest and southbound traffic from Quanzhou. This is to be expected with the northeasterly winds of the winter monsoon. Also, the optimal wind patterns were such that our digital navigator with its generalized ship of the era is able to match the performance depicted on the Selden map. In contrast, many of the simulated routes cut through the open South China Sea instead of following the coasts of China and Vietnam (Fig. 11a). Perhaps the most striking deviation is the route to Brunei, which goes through the Spratly Islands despite the extra cost assigned to their proximity. West of the Mekong, the simulated routes resemble the Selden map more, but this is partly due to geographical limitations.
The route Quanzhou-Mekong Island was chosen for detailed analysis because of its northeast-southwest nature, suitable length for processing and special interest to the "The ports and harbours of Southeast Asia: Human-environment entanglements in Early Modern maritime trade networks" project (Walker Vadillo, 2020). The simulation was repeated daily throughout the year for both directions and the results are presented in Fig. 9. Based on the winds of 1979 and early 1980, the digital navigator could reach Mekong Island in a reasonable amount of time -c. 10 to 30 days -if setting sail from January to mid-April and again from mid-September through December. The shortest simulated duration of 7.7 (7.75) days was found on the launch day of 10th of December and the fastest month on average was November with 9.3 (9.4) days. On the opposite route from Mekong Island to Quanzhou, the window for favourable launch dates was shorter lasting from early May to mid-August. The quickest voyage of 9.9 (10.0) days was started on 9th of August, and June was the quickest month to set sail on with an average duration of 18.4 (19.0) days.
Nevertheless, when considering two-way voyages, the fastest month might not be the optimal one to start the voyage: departing closer to the change in the monsoon winds can save total time even if the one-way voyage is somewhat slower (c.f. e.g. Miksic, 2013: 37-38;Reid, 1988: 64-65). In addition, the time needed for handling the cargo needs to be considered. For example, in the mid-nineteenth century, the sailing between Shangai and Antung took 6-7 days, but the number might be as high as 20 days with loading and unloading included (Yoshino, 1979: 167). If we presume the sailing durations for the following year to be identical, afford the digital navigator a 30-day stay at Mekong Island and sum up the durations, then the launch days from Quanzhou in January to March become indeed more attractive as can be seen from Fig. 10.
During the favourable sailing seasons, the winds were rather dependable. In contrast, the sailing seasons end abruptly: stark jumps can be seen in the sailing durations for the ships launching in mid-April for the Quanzhou-Mekong Island route and in mid-August for its Mekong Island-Quanzhou counterpart. The switches between the winter and summer monsoons happen very quickly in mid-May and mid-September, and setting sail a few days too late could mean considerable delays. Here it is important to note that in adverse winds the routes generated by the LCP analysis can become quite adventurous especially when there is plenty of room to manoeuvre: a small example can be seen in Fig. 11c, where ships launched on 1st-4th of June make a dash north towards the Chinese coast only to loop back when more optimal wind conditions open up elsewhere. Because the digital navigator can consider movement to the eight adjacent cells only, the sector for close-hauled being in this case narrow and the no-go zone wide, and the speed in other points of sail being high in comparison, it has a tendency in headwinds to move sideways or even backtrack in search of more favourable winds rather than simulate tacking properly.   Fig. 9 Daily simulated sailing durations on the routes Quanzhou-Mekong Island (top) and Mekong Island-Quanzhou (bottom). The accumulated cost is shown in darker and loop-based duration in lighter grey (see "Macro 2: Digital Navigator" for further details). The red line denotes the sailing duration calculated from monthly mean winds. Minimum and maximum winds encountered on the route are marked with green and blue lines respectively This is especially true for the simulated routes with very high durations, and they should not be interpreted as actual routes that any real navigator would take but as periods when sailing was not likely to be attempted. In Fig. 11 b and c, we see again how most of the optimized routes between Quanzhou and Mekong Island run through the open sea rather than follow the coast as on the Selden map. To test how much time is saved by this, the voyages starting in November were simulated again for the route Quanzhou-Mekong Island with the digital navigator confined within 150 km from the coast and major islands by introducing a high extra cost for areas further away. The results are shown in Figs. 11d and 12. On average the coastal route was only 0.9 (1.0) days slower and could have been preferable e.g. for the ease of navigation, access to other ports along the route, and/or safety considerations. Also, in reality the gap could be even smaller due to the current flowing along the Chinese and Vietnamese cost during the winter (see "Method").
The macro was programmed to record also the highest and lowest winds encountered during each voyage. If 45 knot windsi.e. mid force 9 -is set as the limit when ships are lost in storm, ships launched from Quanzhou on the 1st of August and the 7th of October would suffer that fate; for example, the latter ship gets caught in the typhoon Sarah. On the return voyage, a single ship would be lost on the 11th of July. Altogether in 1979 there were 10 storms ranging from tropical depressions to typhoons that could affect the traffic between Quanzhou and Mekong Island, all occurring between May and October (Royal Observatory Hong Kong, n.d.: frontispiece). Based on data from 1991-2011, June through November are the most likely months on average for typhoons in the South China Sea (Haghroosta & Ismail, 2017: 30-31, fig. 3). The mariners of the past were well aware of these types of patterns, and it could also influence their decision-making on when was the best time to start a voyage (e.g. Arasaratnam, 1986: 34-37).
In addition, monthly mean winds were calculated from the 6-hourly data and LCP analyses run on them for comparison (Fig. 9). In this particular case, the monthly durations in favourable winds correlate fairly well with the results from the segmented analyses -most likely due to the stableness of the winds -but the correlation breaks down during the unfavourable months. Also, the monthly analysis is "blind" to changes that span over months: e.g. on the route Quanzhou-Mekong the  ships departing in September can utilize the fairer winds in October but this naturally is not reflected in the monthly mean LCPs.

Conclusions
The main aim of this article was the development of the sequential LCP method and the macros needed to run, and it proved successful. The method was tested in the South China Sea region against the Selden Map of China, where the simulated sailing durations matched the values calculated from the map reasonably well. The wind patterns of the area are governed by the biannually changing East Asian Monsoon, and the results demonstrate that the winds are capable of propelling a generalized ship of the era at the pace set by the Selden map during the favourable sailing season. The more detailed daily analyses of the routes Quanzhou-Mekong Island and Mekong Island-Quanzhou show that the monsoon winds guarantee quite dependable sailing during the favourable months as expected. On the analysed year the window for southwesterly travel along the winter monsoon winds was longer than into the opposite direction, and sailing northeast was found to be somewhat slower than to the southwest. The fastest sailing durations on average were found on the Quanzhou-Mekong Island route for the launch dates in November, yet choosing the time to sail was not that straightforward and there were also other considerations in play: avoiding the typhoon season and minimizing the wait for the monsoon winds to change for the return voyage could make launch dates in January to March more attractive options. However, when interpreting the results from the LCP analysis, it is important to recognize that it excels in finding the optimized routes in favourable conditions but cannot in many cases properly simulate tacking in adverse winds. The sequential LCP method provides a new option to the modelling of directed sailpowered seafaring of the past based on deterministic time series wind data. In the South China Sea region, the model can be used to estimate the sailing durations of routes not present on the Selden map or other rutters. It can offer insights into the cyclical nature of seaborne movement and trade in Southeast Asia and how short-term changes in the wind patterns affect the optimal sailing routes and durations. The model can also be refined further if experimental performance data for historical junks becomes available.
With minor modifications to the macro, the method can be applied to similar research questions in other time periods and geographical settings such as the Mediterranean, Caribbean or the Baltic Sea regions. Moreover, as Slayton (2021: pers. comm.) has remarked, comparing models against each other could prove to be a fruitful avenue of research. In this case, a juxtaposition of my results to Davies and Bickler's (2015) tool or to modern navigation software (Gal et al., 2021;Meltemus, 2021) would be extremely interesting: due to their design, the effects of tacking and ocean currents could be tested with them on the Selden map routes. As we move forward with these approaches and begin to fine tune our results, it is possible to get a more accurate reading of travel times by sea, which, combined with further advancements in our knowledge of nautical technology, can help us get a greater picture of how maritime networks actually worked. By gaining a better understanding of how ships moved through the seas, we can begin to think about time, and how this somewhat elusive component affected human interactions in the past. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.