Evaluation of the Heterogeneous Spatial Distribution of Population and Stores/Facilities by Multifractal Analysis

The spatial distribution of a population is not homogeneous—some areas attract many residents, while others do not. The spatial distribution of stores and facilities that have been coevolving with that of the population is also heterogeneous. Previous studies have shown that multifractality is a characteristic of the spatial distribution of a population, as well as other quantities associated with the urban system. We found that stores/facilities belonging to some categories also exhibit multifractality in spatial distribution. We quantified the spatial distributions of the population and stores/facilities in each category by multifractal analysis and compared their multifractal properties. Multifractal measures that reflect the heterogeneity of the densities in each location were able to capture additional features that cannot be seen when only the box-counting dimension was observed. Further, high concentrations of stores/facilities in categories relating to professional or commercial businesses were observed, consistent with previous studies on the scaling law, another pattern observed in urban morphology. We discuss the implications of the multifractal properties on the arrangement of locations where stores/facilities are concentrated. We believe that multifractal analysis is a powerful tool for the quantification of spatial distributions and expect that our interpretation of multifractal measures will stimulate further investigations into urban morphology.


INTRODUCTION
Fractality has been observed in various spatial distributions relating to the morphology of cities, e.g., population, buildings, land price, and street networks [1][2][3][4][5][6][7][8][9][10]. Fractality is represented by the nature where the mass (e.g., the population) in a region exhibits power law dependence on the size of the region. The power-law exponent is called the fractal dimension [11]. The abovementioned spatial distributions also exhibit heterogeneity in fractality in the sense that the locally measured fractal dimension around each spot diverges in the structure, known as multifractality [4,12,13].
The location of stores/facilities should depend on and, in turn, affect the spatial distribution of the population. People may choose the location of their settlement based on the availability and variety of stores/facilities. On the other hand, companies may invest in the construction of stores/facilities to secure employers and customers [14]. It is known that the nature of agglomeration of stores/facilities depends on their type. For example, previous studies showed that facilities relating to businesses offering professional services tend to concentrate in areas with a large population [14][15][16]. One of these studies evaluated the concentration by investigating the scaling law between the population and densities of facilities, showing that the number of facilities in a city increases with a power of the population in the city [16]. The scaling law is a universal pattern observed in urban morphology, e.g., accessibility, road surface, and crime [16][17][18][19][20][21]. It should be interesting to see whether the agglomeration of stores/facilities exhibits multifractality, another urban-related feature.
If the spatial distribution of stores/facilities exhibits multifractality, multifractal analysis can be performed on the structure and is a strong tool to quantify the feature of it. In multifractal analysis, the resulted curve exhibits the correspondence between the locally measured fractal dimension, called singularity strength α, and the fractal dimension of the arrangement of spots that exhibit the singularity strength α, called spectrum f (α) [4,13].
Besides multifractal analysis, various methods have been applied to evaluate spatial distributions including that of the population and stores/facilities. Average nearest neighbor (ANN), which is the averaged distance to the nearest point, is an indicator of the spatial clustering [22][23][24]. The distribution of nearest neighbor distance is also used to evaluate the qualitative feature of spatial structure [24,25]. On the other hand, when grid lines are drawn on the space, probability of finding neighboring cells both of which are occupied is often calculated to assess the degree of clustering [26,27]. Regarding the industrial coagglomeration between two industries, the extent to which the facilities of these industries are in the same region is evaluated [14,28]. Compared to these methods, an advantage of multifractal analysis should be that it can demonstrate both local and global features of the spatial distribution [9]. By multifractal analysis, the strength of local concentration can be captured by the singularity strength. The global view of the arrangement of spots with each level of concentration, on the other hand, is incorporated into the spectrum.
In this study, we aim to (1) determine if the spatial distributions of stores/facilities in various categories exhibit multifractality, and if so, (2) determine the characteristics captured by the multifractal properties of each spatial distribution.
We chose the largest metropolitan area in Japan, the Kantō area, as the object of our analysis. We investigated the multifractal properties of the spatial distributions of the population and stores/facilities in Kantō. Our analysis showed that the spatial distribution of stores/facilities in some categories exhibit multifractality, as does that of the population. Though these spatial distributions are on the same geographical substrate, their multifractal properties are significantly different from each other. This paper is organized as follows. The principles of fractal geometry, multifractality, and the generalized dimension, as well as the methods of our analysis are presented in section 2. We show the results of multifractal analysis of the spatial distributions in section 3. We also compare the multifractal properties of these spatial distributions and highlight the information extractable from the multifractal measures of the spatial distributions. We discuss the results in section 4.

Data
Data for the spatial distributions of the population and stores/facilities were extracted from the Japanese 100-Meter Estimated Mesh Data of the 2015 National Censuses (Zenrin Marketing Solutions Co., Ltd.) and the Corporate Telephone Directory Database Telepoint with Coordinates (Telepoint Pack! provided by ZENRIN Co., Ltd.) of 2017, respectively. The former dataset contains data on the estimated population in each mesh. The length of each side of a mesh is ∼100 m, while the exact size is 3 s in the latitude direction and 4.5 s in the longitude direction. The latter dataset contains the geospatial information of each store/facility. Stores/facilities are categorized hierarchically. In the largest classification, which we adopted, there are 39 categories as shown in Table 1.
The analyzed area is a part of Kantō in Japan, that includes the capital, Tokyo, and a major industrial area, the Keihin industrial area. The range of the latitude is from 35 • 29 ′ 54 ′′ to 35 • 55 ′ 30 ′′ and that of the longitude is from 139 • 16 ′ 52.5 ′′ to 139 • 55 ′ 16.5 ′′ (Figure 1). There are 2 9 × 2 9 (= 262,144) meshes inside this region. Table 1 shows a summary of the data analyzed in this study. The total population (the total number of stores/facilities) in the analyzed area and the number of non-empty meshes are shown in the third and fourth columns, respectively. Here, the meshes/boxes with non-zero populations (stores/facilities) are called support. The maximum and the mean population (number of stores/facilities) in a mesh on the support are shown from the fifth to the sixth columns. Figure 2A shows the spatial distribution of the population in the analyzed area. The logarithm of the proportion of the population in each mesh to the total population is represented by the heatmap. The other panels (B-L) show the spatial distributions of stores/facilities in 11 categories.

Multifractality
We briefly introduce the concepts of fractal geometry and multifractality. There are several ways of defining fractals and multifractals; we present one of them here [9,[11][12][13]29]. Additionally, we limit our explanation to structures embedded in R 2 , while higher dimensions of fractal and multifractal structures can be generally defined.
When a mass (m(ε)) composing a structure within a region of size ε increases with ε according to m(ε) ∼ ε D , the structure is regarded as having fractal characteristics. Here, the "size" is, for example, the length of a side if the region is a square. The exponent D is the fractal dimension of the structure. A more precise definition of the fractal dimension of a structure X is the one by the following box-counting method [11]. Let us assume that a structure X is covered with boxes of size ε. Let N(ε) be the minimum number of such boxes required to cover the structure. The star beside the ID indicates that the spatial distribution of that category exhibits multifractality. The star-marked categories also have their abbreviation enclosed in brackets. The total population (the total number of stores/facilities) in the analyzed area is shown in the third column. The number of meshes that is not empty is shown in the fourth column. These non-empty areas are called support. The maximum and the mean population (number of stores/facilities) in a mesh on the support are shown in the fifth and the sixth columns, respectively.
Then the box-counting dimension D is defined as: We introduce multifractality. We again consider a set X and a function µ on X that gives a quantity, such as the density, at each point x ∈ X. Let us assume that X is divided by square boxes that have the same size ε. For the i-th box of size ε, C i,ε , the value P i,ε is called the probability measure on the box: If P i,ε and ε have the following power-law relationship for any i: then fractality can be seen around each point of X. Here we assume that the exponent α diverges in X and let N(ε, α) be the number of boxes that satisfy the set X can be regarded as having a multifractal structure. Exponent α, which can be regarded as the local fractal dimension, is called the singularity strength. On the other hand, exponent f (α) stands for the box-counting dimension of the set of points with the singularity strength α. This dimension f (α) is called the spectrum. In this paper, the curve of (α, f (α)) is called the multifractal curve. Figure 3 is an example of the relationship between α and f (α) in the spatial distribution of the population (Figure 2A). Each panel in Figure 3 shows the units that exhibit the singularity strength α within the range shown on each panel.
We derived the singularity strength of each unit by estimating the exponent in the following relationship P i,ε ∼ ε α (see Equation 3). The box-counting dimension of the arrangement of units with the singularity strength α, f (α), was also derived based on Equation (4)

Generalized Dimension
We introduce the generalized dimension and explain the relationship between the singularity strength, the spectrum, and the generalized dimension [29,30]. The q-th generalized dimension D q is defined as follows. First, we define τ q as Then the generalized dimension D q for q = 1 is defined as Frontiers in Physics | www.frontiersin.org In the case of q = 1, In the summation on the right-hand side of each of Equations (5), , let the i-th term be summed when the i-th box is not empty, i.e., P i,ε = 0. Here, the generalized dimension is equal to the box-counting dimension of the support D when q is zero.
It is known that the following values of α q and f (α q ) give the approximation of the pair of α and f (α) for each q:

Method of Analysis
In this study, we are concerned with finite morphologies that do not contain infinitesimal structures as the smallest unit of our data is a 100-m mesh. Therefore, we cannot rigorously calculate Here the size of each unit is 2 3 /2 9 on one side, and α is derived based on Equation (3). Also, the spectrum f(α) is derived based on Equation (4), only when the power-law relationship in Equation (4) holds.
the generalized dimension, singularity strength, and spectrum. Instead, we consider the range of q and ε in which the structure can be regarded as exhibiting (multi)fractality, by evaluating the range of q and ε in which the following relationship holds: Thus, we examine the linearity of the relationship between log i P q i,ε and log ε by the frequently used method [3,4,7,8,[31][32][33][34]. Then we regard d as τ q by Equation (5) In our analysis, we assign the size of one side of a 100-m mesh to ε = 1/2 9 . Grid lines are drawn on the analyzed area so that it is covered by non-overlapping boxes with the same size [3,7]. As the size of boxes, we consider ε = 1/2 9 , 2/2 9 , 2 2 /2 9 , ..., 1. We define the probability measure P i,ε on the i-th box with the size ε as the proportion of the population (the number of stores/facilities) in the i-th box to the total population (the total number of stores/facilities).
Subsequently, for the population, we evaluate the linearity of the relationship between log i P q i,ε and log ε for various ranges of ε and q as shown in Figure 4. Regarding the criterion for this linearity, we examined whether or not the coefficient of determination of the linear regression, R 2 , exceeds 0.99. A linear relationship was not observed when q takes a negative value and when the range of ε includes values < 2 3 /2 9 . Consequently, we considered the range of ε from 2 3 /2 9 to 2 9 /2 9 , and the range of q from 0 to 20. We performed a multifractal analysis on these ranges of ε and q for the spatial distribution of the population. Furthermore, for the spatial distribution of stores/facilities in each category, we examined multifractality on the same ranges of ε and q as that for the population. Plots of log i P q i,ε against log ε for all categories of stores/facilities are shown in Figures S3, S4. The star marker was added beside the ID in the first column of Table 1 if the category showed multifractality in the spatial distribution. Also, the spatial distribution of such categories are shown in Figures 2B-L. In this figure and in the following discussion, abbreviations are used for these categories-the abbreviation of each category is enclosed by brackets after the name of the category in the second column of Table 1. We obtained τ q as the slope of the linear regression of log i P q i,ε by log ε, and derived D q from Equations (6, 7). To derive the singularity strength α q and the spectrum f (α q ), we used the following formulae: where µ i,ε,q = P q i,ε / j P q j,ε . These formulae can be directly derived from Equation (5), (8), and (9) [12,35,36]. We obtained α q by the linear regression between i µ i,ε,q log P i,ε and log ε, and obtained f (α q ) by the linear regression between i µ i,ε,q log µ i,ε,q and log ε.

Density of Units and Multifractality
As we briefly mentioned in section 2.4, multifractality could not be observed for the spatial distribution of the population when the value of q is negative. Knowledge of the relationship between the densities in a unit, the value of q, and the multifractal measures is helpful for interpreting the results of the multifractal analysis. In Equation (5) for τ q , a larger value of q corresponds to a greater contribution of the boxes with large probability measures P i,ε to the sum. Therefore, the boxes with high densities are significantly incorporated into the calculation of the q-th generalized dimension D q when the value of q is large, see Equation (6). Note that no differences in density is considered when q = 0. The pair of the singularity strength and the spectrum also has a relationship with q by Equations (8), (9). Therefore, the multifractal measures of D q , α q , and f (α q ) for a large q, significantly reflect those units with high densities. The fact that multifractality was only observed in the spatial distribution of the population with positive q values, indicates that multifractality cannot be observed in sparsely populated areas. Figure 5 shows the q-th generalized dimension vs. q for the population and stores/facilities in four categories: 20FIN, 25WH, 31AED, and 33HW. The generalized dimension for stores/facilities in the other categories is shown in Figure S5. As mentioned before, the generalized dimension for q = 0, D 0 , is equal to the box-counting dimension of the support. For each category, the value of D 0 is close to 2, which is the dimension of R 2 to which each spatial distribution is embedded. The rate of decline in the generalized dimension with q significantly varies between the categories. The generalized dimension for 20FIN drops dramatically with q, while that of the population decreases minimally with q. In the cases of 25WH and 33HW, the decline of the generalized dimension is milder than that of 20FIN and 31AED. Considering the relationship between the value of q and the densities in the boxes, | D 20 | : = |D 20 − D 0 | can be one of the indicators for the strength of the heterogeneity in spatial density distribution, and 20FIN and 31AED can be regarded as having strong heterogeneity in spatial distributions.

Multifractal Curve
The multifractal curves are shown in Figure 6A. Before we investigate each multifractal curve, let us revisit the relationship between the densities in a unit and the value of q, and the singularity strength α q [9]. Recall that each point on the multifractal curve can be derived for each q (Equations 11, 12). Generally speaking, in each multifractal curve, the greater is the value of q, the lower is the value of α q , i.e., the plot of a multifractal curve proceeds to the left as the value of q increases. Therefore, by the relationship between the value of q and the densities in the units, a unit with a high density tends to have a small singularity strength-we can observe this by comparing Figure 2A with Figure 3. An interpretation is that around a unit with a density significantly higher than the surroundings, the probability measure P i,ε does not increase rapidly by expanding the area, and vice versa (see Equation 3). Therefore, in Figure 6A, multifractal curves with a significantly low singularity strength α q , such as 20FIN and 31AED, indicate the existence of units with a density much higher than the surroundings. Also, the spectrum f (α q ) is quite low for small singularity strengths α q in the cases of 20FIN and 31AED. Recall that the spectrum f (α q ) represents the box-counting dimension of the arrangement of units with the singularity strength α q . The pair of small α q and f (α q ), thus, is an indication that there are a few isolated units, each with an extremely high density.
Multifractal curves contain information not only at small singularity strengths but also across the entire range. The multifractal curve of 13MP has similar α q and f (α q ) values with that of 25WH, 29SF, 35LS, and 37SE when α q is quite small. However, in the mid-range of α q , the value of f (α q ) for 13MP is significantly greater than that of the others. Therefore, in the case of 13MP, units with a mid-range density (compared to the surroundings) remain concentrated compared to the case of the other categories. This may correspond to the yellowish units gathered on the right of panel (C) in Figure 2. On the other hand, a feature of the multifractal curve of the FIGURE 6 | (A) Multifractal curve. The pair of the singularity strength (abscissa) and the spectrum (ordinate), (α q , f(α q )), for each q of the population and stores/facilities. Each multifractal curve is also separately shown in Figure S6.
population is that neither α q nor f (α q ) declines rapidly with q. This corresponds to the following features of the spatial distribution of the population: The densities in the units does not widely diverge and the arrangement of the units with each density does not change dramatically with the density (Figure 2A). Figure 6B shows the plots of f (α q ) : = f (α q ) − f (α 0 ) against α q . Recall that f (α 0 ) is the box-counting dimension of the support. In the value of f (α 0 ), the differences in the densities in the units is not incorporated. Thus, the vertical axis f (α q ) represents the degree of decline in the spectrum from the box-counting dimension of the support when the differences in the densities in the units is gradually emphasized. Another intuitive meaning of f (α q ) is the nature of the difference in arrangement between the panels in Figure 3. In the case of 31AED, in which the pair of a small α q and f (α q ) exists, we observe that the value of f (α q ) is also small for small α q . This indicates that the arrangement of units with small α q is quite sparse compared to the support. Therefore, we can infer the following characteristics of the spatial distribution of 31AED: There are a few isolated locations where stores/facilities in 31AED are extremely concentrated, while covering a vast region. Figure 6C shows plots of f (α q ) against α q : = α q − α 0 . The horizontal axis shows that how the singularity strength declines as q increases from the one calculated under the situation where differences in the densities in the units were not taken into account. The decline of both α q and f (α q ) with q of 29SF is milder than the other categories excepting the population, suggesting that the spatial distribution of stores/facilities is homogeneous in 29SF. Specifically, by comparing with Figure 6A, the value of α q is small across the range of q including q = 0, in the case of 29SF. This may represent the sparsely but relatively uniform scattering of stores/facilities in 29SF ( Figure 2G). Additionally, though the difference in the multifractal curves of 21RE and 32RES is ambiguous for large q in Figure 6A, we can observe that the singularity strength of 21RE has a wider range than 32RES (Figure 6C). When we restrict our observation to the concentrated units, the arrangement of the units in 21RE has a similar box-counting dimension to that of 32RES. On the other hand, in 21RE, we should be able to observe local regions where the density increases mildly with the size of the region, that may not be observed in 32RES.
Finally, Figure 8 shows a summary of the multifractal properties of the spatial distributions of the population and stores/facilities for all categories that have multifractality in the ranges of ε and q tested in this study. In both panels, the horizontal axis shows f (α 0 ), i.e., the box-counting dimension of the support. In f (α 0 ), the differences in the densities in the units is not incorporated-only whether or not each unit is empty is taken into account. The box-counting dimension of the support, of course, can capture a feature of each spatial distribution. The value of f (α 0 ) for each category is near two, that is the dimension of R 2 , but it diverges a little. We can see that f (α 0 ) tends to be small when the spatial distribution is sparse, e.g., 20FIN and 29SF, in (Figures 2D,G, 8).
On the other hand, the vertical axes in Figure 8 show the values that incorporate the differences in the densities in the units. In Figure 8A, the vertical axis shows the difference between the largest and the smallest singularity strength | α 20 |. The color of each marker represents the value of | D 20 |, for each category. The brighter marker color corresponds to the larger value of | D 20 |. In Figure 8B, the vertical axis shows the difference between the largest and the smallest spectrum | f (α 20 )|. The color of each marker again represents the value of D 20 .
The values of | α 20 | and | f (α 20 )| for some categories can diverge even when f (α 0 ) takes almost the same value. For example, | α 20 | of 29SF is much smaller than that of 20FIN, while f (α 0 ) of both of these categories are around 1.86. This result indicates that both 29SF and 20FIN have sparse spatial distributions, but the heterogeneity of the densities in the units for 20FIN is stronger than that of 29SFthis can be observed in Figures 2D,G. In Figure 8A, the plots of 3CC and 33HW are nearer that of the population than the others. This suggests that the nature of the spatial distribution in these categories is similar to that of the population when we consider the large region covered by stores/facilities and the small differences in the relative densities in the units.

DISCUSSION
In this study, we evaluated whether the spatial distributions of the population and stores/facilities exhibit multifractality. Multifractality in the spatial distribution of the population has been demonstrated in previous studies [3,5,10]; we also demonstrated this result in the Kantō area, Japan. However, we were not able to observe multifractality in the population distribution for negative values of q. In the previous studies that evaluated multifractality in city morphologies, the authors carefully examined the range of q in which multifractality was observed [4,7,8]. These studies showed that multifractality can be observed at both positive and negative values of q. As we mentioned above, positive (the negative) values of q corresponds to boxes with high (low) densities. Therefore, previous studies observed multifractality in both densely and sparsely distributed regions. However, they also showed asymmetry in the positive and negative ranges of q and discussed the structural differences in dense/sparse regions. In our case, the range of q for multifractality suggests that the sparsely populated region does not have a structure characterized by multifractals, attributable to the geographical characteristics of the examined region in this study. The examined region includes the mountainous areas in the upper and left sides of each panel in Figure 2 as it represents the general feature of the Japanese terrain. It also contains Tokyo bay in the lower-right corner of each panel (Figure 2). Along Tokyo bay, there are rich residential areas as well as plenty of facilities in various industries and numerous stores. This examined region should be an interesting object to investigate considering these geographical features. However, complex substrates may restrict the range in which multifractality appears.
We also investigated which category of stores/facilities shows multifractality for the same ranges of q and ε as that for the population. Stores/facilities in some categories also exhibit multifractality in the spatial distribution, but the determined multifractal measures significantly depend on the category. Diverging multifractal properties can reflect qualitative differences in the spatial distributions-stores/facilities are sparsely and uniformly scattered in some categories, while others are centralized, e.g., the stations in some categories. Importantly, our analysis showed that the box-counting dimension performs poorly in capturing qualitative diversities in the spatial distributions of these categories. The box-counting dimension captures the arrangement of the units that are not empty. On the other hand, multifractal measures can represent the arrangement of units with a certain density. Multifractal curves can indicate, for example, the existence of units with comparably high densities by evaluating the range of the singularity strength α, and the spatial distribution of such units by the spectrum f (α).
The spatial distribution of the population can be characterized by the high box-counting dimension of the support and the homogeneity measured by the generalized dimension and the multifractal curve. The population is distributed across all the regions examined in this study, while the densities in a unit do not vary significantly. In addition, we will observe a similar spatial distribution, even when we change the filter on units according to these densities. Considering its high box-counting dimension and the homogeneity seen in the range of the singularity strength, the spatial distribution of stores/facilities in 33HW exhibits similar features with that of the population. These stores/facilities cover a large area and the heterogeneity of the densities in the units is low. On the other hand, in the cases of 20FIN and 31AED, the singularity strength and the spectrum sharply decline with q. This shape of the multifractal curves indicates a strong centralization of stores/facilities to a few locations. Some of these concentrated locations presumably correspond to large stations in the capital.
In addition to multifractality of the spatial distribution, the scaling law is also a universal pattern observed in city morphology [16-21, 37, 38]. For example, the population X and a quantity Y related to the city morphology have the relationship of Y ∼ X β . The large scaling exponent β in this relationship, i.e., the super-linear increase of Y with X, indicates a strong concentration of these urban-related objects to locations with a large population. Previous studies also found the scaling law in (industrial) agglomerations and showed that facilities, outputs, and jobs concentrated stronger in cities with a large population when these objects are in a category associated with professional and complex skills or with commercial facilities than when the type is relatively primitive or public-related [16,[19][20][21]. Many of our results are consistent with these previous studies. Our multifractal analysis indicates the centralization of facilities in 20FIN, which is a category related to services requiring professional skills and frequent communication with customers [14]. Also, 31AED, a category related to commercial activities, showed a strong concentration. The concentration of objects related to construction and healthcare was shown to be mild in previous studies, which is also consistent with the mild decline of the singularity strength with q in our results (3CC and 33HW).
While consistency between multifractality and the scaling law exists, an advantage of multifractal analysis should be the richness of information in the result. As we discussed so far, we can quantify the nature of the divergence of the densities in each location and the spatial distribution of each density, by multifractal analysis. For example, for the centralization of stores/facilities in 20FIN and 31AED, multifractal properties further explain the following difference. Considering the large range in the singularity strength and the high box-counting dimension, we can expect to see stores/facilities in 31AED everywhere, with some centralized locations. On the other hand, facilities in 20FIN are encountered only in concentrated locations, which is represented by the overall small singularity strength. Therefore, the results in our analysis exhibit not only the existence of concentrated areas but also the various state of concentration. In this study, we attempted to interpret the multifractal curves that correspond to qualitatively diverging spatial distributions. We hope that our discussion will contribute to future investigations on spatial distributions by multifractal analysis. Additionally, the characteristic of concentration of stores/facilities in each category, which was revealed in this study, should be considered in the actual urban design. For example, the stores/facilities in 31AED is expected to have a tendency to concentrate strongly. Such a tendency of agglomeration should be taken into account in advance when it is required to avoid an extreme concentration of buildings in a landscape.
The temporal development of various urban morphologies, e.g., the spatial distributions of streets and buildings, have been discussed in previous studies [2,4,[7][8][9]. Some of them revealed that the spatial distribution was developed to the packed state and to exhibit features close to monofractals [2,4]. We are also interested in how these developments depend on the category in which the stores/facilities belong. As a future perspective, the comparison of such developments is possible by quantification with multifractal analysis. Furthermore, we expect that the classification of cities is possible by comparing such developments between cities.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because these are available under the contract. Requests to access these datasets should be directed to Zenrin Co., Ltd., https:// www.zenrin.co.jp/product/category/gis/contents/teldt/index. html for Corporate Telephone Directory Database Telepoint with Coordinates (Telepoint Pack!) and to Zenrin Marketing Solutions Co., Ltd., https://www.zenrin-ms.co.jp/gis_marketing/ database/statistics/census_100m_mesh/ for Japanese 100-Meter Estimated Mesh Data of the 2015 National Censuses.