Evaluating the use of publicly available remotely sensed land cover data for areal interpolation

Abstract

Areal interpolation is used to transfer attribute data between geographically incongruous zonal systems. Remotely sensed land cover data are widely used in intelligent areal interpolation methods to solve this problem. This article examines the usefulness of different publicly available remotely sensed land cover data sets as ancillary data used in conjunction with different areal interpolation methods. Two land cover data sets were compiled at the national scale; one by the Multi-Resolution Land Characteristics Consortium (the National Land Cover Dataset or NLCD) and one by the Coastal Change Analysis Program. A third land cover data set was compiled at a regional scale for the state of Connecticut by the Center for Land Use Education and Research. Results show that for areal interpolation, greater detail in the classification of developed areas was important whether the data were developed for use at a national or a regional scale. Even more important is the further enhancement of remotely sensed land use categories by incorporating local road or parcel data layers. The worst performing interpolation method using enhanced remote sensing-derived land cover data produced more accurate results than the best performing method using only the original land cover data. The results also show that parcels produce better enhancements than road buffers because they remove the areas of the roads themselves from population consideration.

Keywords:

• areal interpolation
• remote sensing
• land cover
• NLCD
• C-CAP

1. Introduction

Areal interpolation techniques have been developed to enable data collected for one areal delineation of geographic space to be transferred to a different delineation of geographic space (Lam 1983). This is a type of “change of support” problem, in which the set of regions with known attribute values are identified as source zones, and the set of regions needing estimated values are target zones (Goodchild and Lam 1980). Fotheringham and Rogerson (1993) categorized and discussed four general issues and problems associated with spatial analysis: the modifiable areal unit problem (MAUP), the boundary problem, the spatial sampling problem, and the spatial autocorrelation problem. Areal interpolation is especially associated with the MAUP, which includes the sensitivity of spatial analysis results to the zoning system used to aggregate spatial data (Openshaw and Taylor 1981). Many spatial data sets are aggregated into zones because of privacy issues and storage efficiency, but the diversity of applications, variations in geographic scale, and the spatial distributions of the phenomena under investigation have worked against the unification of areal units into a single standardized system (Visvalingam 1991). Areal interpolation can be used to estimate missing data layers, transfer data from one geography to another, and estimate attribute values for new polygonal units resulting from GIS overlay operations (Mrozinski and Cromley 1999). Thus, areal interpolation can be used to rectify disparate spatially aggregated data layers into a common scale with common boundary delineations before proceeding with further spatial analysis. In addition, the areal interpolation of population is needed in many other applications (Qiu, Zhang, and Zhou 2012).

Overtime, different methods have been developed to improve the efficiency and accuracy of areal interpolators. Different ways for classifying these procedures can be found in the literature, based on differences in model assumptions and data types. Each method makes an assumption regarding the distribution of data within the spatial units that can result in different levels of interpolation accuracy. Methods that use ancillary data have been found to improve the results of the interpolation (Langford 2007). This is because the ancillary data used are often spatially correlated with the phenomenon under investigation in the source zone. Having a better understanding of the distribution of the data in the source zone can improve the inherent disaggregation of the data during the areal interpolation process. The most widely used classification scheme divides areal interpolation methods into two categories: simple methods, which assume data are uniformly distributed within source zones, and intelligent areal interpolation methods, which use ancillary data to infer a spatial distribution of the data of interest (Okabe and Sadahiro 1997; Langford, Maguire, and Unwin 1991; Langford 2006). Intelligent methods are further subdivided into those based on dasymetric mapping (Wright 1936) and those based on a statistical method such as regression.

Another classification is based on the underlying data model used during the interpolation procedure. For a vector-based model, the spatial data used for areal interpolation, including source zones, target zones, ancillary zones, if necessary, and all other spatial zoning systems used in the intermediary steps, are all vector data. However, if raster data are included in any step in the interpolation procedure, the method can be considered as raster-based, even though the source layer is usually in a vector format. If the target layer is in a vector format, a raster-based model would use vector-to-raster and raster-to-vector operations as intermediary steps and the interpolation itself is calculated over a grid. The values for the target zones are simply the aggregation of grid values within each areal unit. If the target layer is a raster, then only an initial vector-to-raster conversion of the source layer is necessary as an intermediate step and the interpolation is performed over a grid. Using remotely sensed data as an ancillary layer in the interpolation process also necessitates certain vector/raster conversions.

The purpose of this analysis is to investigate how useful publicly available remotely sensed land cover data sets are with respect to the process of areal interpolation. The next section reviews different areal interpolation procedures including those used in the analysis.

2. Areal interpolation methods

The original vector-based areal interpolation method is areal weighting, which is based on the geometric overlay of the source and target zones (Goodchild and Lam 1980; Lam 1983). It assumes a uniform density of data values within each source zone. Values are estimated as averages that are proportionally weighted by the ratio of the area of each intersection polygon to the area of the source polygon that contains it. Areal weighting is the most widely used method among all other methods because of its intuitively simple theory, low data and computation requirements, and easily implementation in GIS software (Xie 1995). This method is not used in this study because it does not use any ancillary data.

The binary dasymetric method (BDAS) (Fisher and Langford 1996) is a direct extension of areal weighting in which an ancillary data layer containing a single-valued correlate variable is used. BDAS can also be easily implemented in GIS, which can readily integrate various spatial ancillary data (Wu, Qiu, and Wang 2005). In this method, the area of the correlated variable within the source zones is used, rather than the total area of the source zone, to determine the weighted proportions. A spatially correlated variable that is often used is land cover. In a land cover data set derived from remotely sensed images, a single category such as “developed” is used. All land cover categories that might be associated with population would be compressed into a single category. Target zone population is estimated as:

(1)

where  = the estimated population for the target zone; R_ts  = the area of the correlated variable in the intersection zone between the target and the source zone; P_s  = the observed population in the source zone; R_s  = the area of the correlated variable in the source zone.

A slightly different method is intelligent dasymetric mapping (IDM) (Mennis and Hultgren 2006), which allows multiple classes in the ancillary layer. Although this method uses land cover data, it first converts the data to a vector format before the interpolation step. The method then overlays the source layer and the ancillary layer and the resulting intersection zones are used as the target zone. Only one ancillary class is associated with each target zone. The estimated population for the target zone is:

(2)

where  = the estimated population for a target zone; y_s  = the population of a source zone, which contains the target zone; A_t  = the area of the target zone;  = the estimated density of each ancillary class for the target zone.

IDM incorporates the estimated density of each ancillary class in calculating the proportion of counts in a target zone from that in the source zone. The value of can be specified by the analyst who has some prior knowledge of the density of the ancillary class. Noting that this specification is subjective on the part of the analyst, three objective methods are further proposed to calculate : containment, centroid, and percent cover. All source zones that are completely contained within each individual ancillary class are selected as a sample to calculate the density for that class in the containment method. In the centroid method (IDMC), source zones that have their centroid within each individual ancillary class are selected as the sample to calculate the density for that class. For the percentage method, source zones with areas occupied by a single ancillary class equal or exceeding a threshold percentage, subjectively assigned by the analyst, are selected to calculate the density for that class. Calculating density for each ancillary class for all three sampling methods is:

(3)

where A_s  = the area of one selected source zone; m = the number of selected source zones associated with the ancillary class.

Noting that the densities for each of the ancillary data are globally derived by the ratio of the sum of the counts for all selected source zones to the sum of the areas for all those zones, the authors further described a more local model in which the study area is divided into several regions to account for spatial variation in the relationship between density and ancillary class, with density calculated separately in each individual region. The possibility exists that no source zone is selected for a particular ancillary class. In the containment method, for example, no source zone was fully contained within one particular ancillary class; thus, the density for this class cannot be calculated using the method described above. In this case, they proposed first calculating the raw estimated count for the target zone containing an ancillary class with unknown density using the previously estimated densities of other ancillary classes that are in the same source zone as this particular class:

(4)

where  = the raw estimated count of the target zone containing an ancillary class with unknown density;  = the estimated density of other ancillary classes;  = the area of the target zones containing other classes;  = the area of the target zones containing the ancillary classes with unknown density.

Note that is not the final estimated count for the target zone, rather it is used to calculate the density of the ancillary class with unknown density as follows:

(5)

where  = the estimated density for the ancillary class with previously unknown density; p = the number of target zones in the entire study area containing the ancillary class.

Then, the estimated density for this particular ancillary class is used along with Equation (2)(2) to calculate the final count for those target zones associated with this particular ancillary class. Other dasymetric methods using multiple classes have been developed by Eicher and Brewer (2001) and Langford (2006), but are not used in this analysis.

In the statistical approach to areal interpolation, different regression models have been applied to estimate population density values for each class of the correlated variables. Langford, Maguire, and Unwin (1991) used ordinary least squares (OLS) regression with population count as the dependent variable and pixel counts of individual land cover classes as independent variables to model the relationship between population density and land cover classes. Noting the possibility of negative counts using OLS, Flowerdew and Green (1989) suggested that Poisson regression is theoretically preferable for modeling counts, as negative population estimation could be avoided. They proposed a Poisson model to interpolate count variables using ancillary binary data for target zones. Furthermore, with the existing target zone information, Flowerdew and Green (1991) improved their Poisson model by synthesizing the expectation and maximum likelihood (EM) algorithm, which was originally developed by Dempster, Laird, and Rubin (1977) to solve problems of missing data. Flowerdew and Green (1994) then extended this EM method to handle continuous variables having a normal distribution.

However, the early regression models generally underperformed in comparison with a simple binary dasymetric model (Fisher and Langford 1995). Population density varied more among different places than among different ancillary classes. To make the relationship between population and an ancillary class in a regression model not only vary over ancillary classes but also over space, Yuan, Smith, and Limp (1997) developed a regional model that regressed population on land cover types in each county of their study area. However, Langford (2006) argued that their method for dividing the study area into sub-regions to fit the local regression model is not based on the distribution of population. The counties used by Yuan, Smith, and Limp (1997) are administrative areas where there is little basis for determining the underlying distribution of population. Furthermore, the variation of local model parameters among counties indicates that some degree of spatial variation exists in the relationship between population density and land cover type. Langford (2006) inferred that the spatial variation in the relationship between population density and land cover type is likely to be continuous, because variation is likely within regions if it exists between regions.

To overcome this problem, Lin, Cromley, and Zhang (2011) proposed using geographically weighted regression (GWR) (Charlton and Fotheringham 2009) for modeling the spatial variation in the relationship between population density and land cover type. Different GWR models were proposed to investigate the influence of bandwidths and the number of classes in the estimation of population density. Because there is usually a residual associated with each source observation when applying GWR (or OLS regression) to areal interpolation, so the population estimates are not volume-preserving. Thus, a scaling step is required to preserve the original population size in the source zone. With scaling, the GWR interpolator outperformed areal weighting, binary dasymetric, and a corresponding global regression interpolator. Also fewer, more relevant classes were found to be more accurate than including a large number of categories.

Another local statistical interpolation model using quantile regression (QR) (Koenker and Bassett 1978) was developed recently by Cromley, Hanink, and Bentley (2012). Theoretically, two advantages were identified for a QR-based interpolator compared with other regression interpolators. One is that every observation can perfectly fit the local regression hyperplane associated with the observation's quantile level; thus, QR-based areal interpolation is inherently volume-preserving and does not require a scaling step. The other is that the estimated coefficients found by QR can easily be constrained to be non-negative and the intercept term can be set to zero; thus, the estimated coefficients can be directly interpolated as a population density for each ancillary class at a specific location. Although the authors did not compare it with GWR, they found that the QR interpolator outperformed areal weighting, binary dasymetric, OLS regression, and a spatial error model-based interpolator (Anselin 2003).

Most of the existing literature compares different areal interpolation methods using a single ancillary data layer such as land cover, parcel data, a road network, or other physical or socio-economic characteristics. Fisher and Langford (2006) also examined the effect of the errors in a classified Landsat image on the areal interpolation of population counts. In their research, the classification errors were simulated by assigning any pixel in a classified satellite image to an alternative random cover type between 0% and 100% at 10% increments. The simulation was further divided into constrained and unconstrained, mainly based on whether the original proportions of the different land cover types were preserved during the simulating process. They concluded that dasymetric method, which is their main test method, is robust against classification errors in satellite imagery.

Instead of testing the simulated errors in classified satellite imagery, this article focuses more on examining the utility of different publicly available preclassified remotely sensed land cover data sets with respect to different areal interpolation methods. Until recently, most of the publicly available preclassified land cover data sets provided only an Anderson Level I-like classification, which distinguishes only among the broadest land cover types: developed agriculture, forest, water, wetland, etc. (Anderson et al. 1976). However, the Anderson Level II classification further divides developed areas into residential, commercial, and industrial areas. Falcone and Homer (2012) noted that derivation of the Anderson Level II classification, based on the spectral information in remotely sensed imagery, can be enhanced by using additional socio-economic data layers. In their research, they further classified the developed pixels in the National Land Cover Database 2006 (NLCD 2006) into very low density residential, low density residential, medium density residential, high density residential, commercial-industrial with low-medium density, commercial-industrial with high density, institutions, transportation, and open space. The main enhancement data they used in their research is the Homeland Security Infrastructure Program (HSIP) HSIP-Gold 2010 data set (http://www.hifldwg.org/), which provides the locations and shapes of institutions, transportation networks, recreation, commercial, and industrial places. The second enhancement data are vector data obtained from ESRI (http://www.esri.com/) containing primarily institutional and recreational features. The third enhancement data set is a historical land-use data set from the Geographic Information Retrieval and Analysis System (Price et al. 2006) and NLCD1992 (Vogelmann et al. 2001). Both data sets delineate several Level II land use classes.

Similar data layers have been used in areal interpolation as alternatives to remotely sensed data. Road buffers (Mrozinski and Cromley 1999), road length (Xie 1995), and parcel data (Maantay, Maroko, and Herrmann 2007; Tapp 2010) have all been used as ancillary data. Other data layers are rarely used in conjunction with remotely sensed data to improve interpolation accuracy. This research not only compares the publicly available preclassified land cover data sets for different areal interpolation methods, but it also examines different enhancements of the land cover data sets for improving the accuracy of these areal interpolation methods.

3. Study area and data

The study area consists of a nine-town region in Hartford County, Connecticut (Figure 1). In the 2010 census, the nine-town study area had a total population of 396,435. These towns have been chosen as the study area because this region is characterized by various types of land use, including residential areas with different housing densities and socio-economic structures, commercial and industrial areas, and open space. In addition, this region contains a variety of population densities, from dense urban centers such as Hartford to sparsely populated rural forested areas. Census tract and block group boundaries data, for the source zones and target zones, respectively, were downloaded from the US Census Bureau's 2010 Topologically Integrated Geographic Encoding and Reference (TIGER/Line) Shapefiles Main Page (http://www.census.gov/geo/www/tiger/tgrshp2010/tgrshp2010.html). Figure 1 shows the study area overlaid with the target and source zones. The population data were extracted from 2010 Census Summary File 1 (SF1) Table P1, and downloaded through the US Census Bureau's web data retrieval interface, American FactFinder (http://factfinder2.census.gov). The geography boundary files and the demographic data were then joined in ArcGIS 10.1.

Figure 1. The nine-town region overlaid with the source and target zones. The nine towns are Bloomfield, East Hartford, East Windsor, Hartford, Manchester, South Windsor, West Hartford, Windsor, and Windsor Locks. One tract in Winsor Locks was not included in the study region because it is Bradley International Airport.

The first public land use/land cover (LULC) data set used as thematic ancillary data was compiled nationally and was based primarily on the unsupervised classification of Landsat 7 Enhanced Thematic Mapper+ (ETM+) satellite imagery acquired in 2006 with a 30 × 30 meters spatial resolution (Fry et al. 2011). This data set is available from the NLCD 2006 from the Multi-Resolution Land Characteristics Consortium (MRLC) (http://www.mrlc.gov/). The NLCD 2006 was created on a Landsat WRS Path-Row basis and mosaicked to create a seamless national product. The land cover classification was performed on an individual path/row basis to make it more local (Xian, Homer, and Fry 2009). The study area is located in Path 13 Row 31, which includes most of Connecticut. The second public LULC data set, which was also derived at the national level, is primarily based on Landsat 5 Thematic Mapper scenes and produced by the Coastal Change Analysis Program (C-CAP). The data set was downloaded from the National Oceanic and Atmospheric Administration (NOAA) coastal service center (http://www.csc.noaa.gov). The third LULC data set, which was derived for use at the state level, is mainly based on Landsat Thematic Mapper satellite imagery and was produced by the Center for Land Use Education and Research (CLEAR) at the University of Connecticut. The data were acquired from CLEAR's website (http://clear.uconn.edu).

Since the satellite imagery, providing the foundation for land cover classification for the three preclassified land cover data sets, was acquired circa 2006, there is a temporal discrepancy between the thematic LULC ancillary data (dated 2006) and the population data (dated 2010). Although there is a discrepancy between the three LULC data sets and the population data, this discrepancy is the same for each LULC data set. All three preclassified LULC data sets were clipped to the study area and registered to the Connecticut State Plane Coordinate System (NAD83) at a 30 × 30 meters spatial resolution.

It is apparent that most of the population is concentrated in the developed area and that some categories of the land cover in the original LULC data set may not be clearly related to the population distribution. As has been noted, having more categories does not necessarily lead to improved results (Lin, Cromley, and Zhang 2011). Only categories related to developed land in the three LULC data sets were used as ancillary data in the analysis. There are three different developed types in NLCD 2006 and C-CAP data sets: low-, medium-, and high-intensity developed. However, there is only a single developed category in the CLEAR data set.

The road networks were downloaded from the US Census Bureau's 2010 TIGER/Line Shapefiles Main Page (http://www.census.gov/cgi-bin/geo/shapefiles2010/main) and clipped to the study area. Finally, parcel data were obtained from the Connecticut Capitol Region Council of Governments (CRCOG). This data layer was first clipped to the study area, and then single-family features, two-family features, and multiple-family features were extracted as residential parcels. The CRCOG parcel data contains 80,022 single-family parcels, 6036 two-family parcels, and 8177 multiple-family parcels in the study area.

4. Research design

The basic problem is to use population counts given for source zones at the census tract level to estimate population counts at the larger scale of block groups. There are broadly four types of areal interpolation methods included in this study: BDAS, IDM, GWR, and QR. Two IDM sub-models with different population density calculation procedures were tested: one calculated population density for each land cover type based on the source zones that had their centroids located in that land use area (IDMC) and the other calculated population density for each land use type based on the source zones where land use appears most frequently (IDMM). The GWR interpolator also had two sub-models with different methods for the scaling procedure. The first scaled the population calculated from the raw estimated coefficient to preserve source zone counts, which is termed GWR with scaled population (GWRSP). The other made all estimated coefficients non-negative first and then scaled the estimated population to preserve source counts; this is termed GWR with scaled coefficients (GWRSC). The QR model had three sub-models. The first two initially permitted negative coefficients as in GWR; the original coefficients were then scaled in a manner similar to GWRSP (QRSP), and the other used the scaling procedure similar to GWRSC (QRSC). The third followed the procedure outline by Cromley, Hanink, and Bentley (2012) by solving the QR with non-negative coefficients in a linear programming model without the need for a scaling step. This is termed the QR with non-negative coefficients (QRNC).

The global accuracy of each interpolation method was evaluated using the root mean square error (RMSE), the adjusted root mean square error (adj-RMSE), the mean absolute error (MAE), and the adjusted absolute error (adj-MAE). They were calculated as follows:

(6)

(7)

(8)

(9)

Where  = the interpolated population for the ith target zone;  = the actual population for the ith target zone; n = the number of observations in the target zones.

RMSE is based on the variance of the estimator; RMSE will increase as the variance of error magnitudes increases (Res, Willmott, and Matsuura 2005). Since the errors are squared before they are averaged, the RMSE assigns relatively high weights to large errors and relatively low weights to small errors. Therefore, the adj-RMSE measure was developed (Gregory 2000) to adjust error magnitudes by the size of the population counts. However, Qiu, Woller, and Briggs (2003) argued that the scaling procedure by population counts for each observation would be more likely to exaggerate the errors for the observations with small population counts and depreciate the errors for the observations with large population counts. MAE is based on the average absolute value of the errors. Res, Willmott, and Matsuura (2005) preferred MAE over RMSE in describing average model-performance error. The adjusted MAE score was also included to reduce the impact of larger population counts.

In the first stage of the analysis, the interpolations used only the extracted land cover categories as ancillary data. Population was assumed to be located within these land cover areas. In the second stage, the land cover categories were enhanced before being used in the interpolation. One enhancement was associated with a 100 feet buffer of only secondary roads (MTFCC S1200), local neighborhood roads, rural roads, and city streets (MTFCC S1400). These roads were selected from the original road network, because certain features such as primary highways, their ramps, and parking lot roads are not likely to have population living alongside. All populations are assumed to reside only within the extracted land cover pixels that fall within this road buffer. The buffer enhancement was performed for all the three original LULC data sets. The second enhancement intersected the single-family, two-family, and multiple-family parcel layer with the three LULC data sets, respectively, to select only those developed pixels within each type of residential parcels.

5. Results and discussion

Table 1 provides the global accuracy of each areal interpolation method using the three different LULC data sets without any enhancement. There is only one developed type in the CLEAR data set, and as shown by Cromley, Hanink, and Bentley (2012), any interpolator based on regression with only one independent variable and no intercept with a scaling step will produce results identical to the binary dasymetric interpolator. All RMSE, adj-RMSE, MAE, and adj-MAE values for GWR- and QR-based interpolators are the same as those for BDAS interpolator under the CLEAR category. IDMC and IDMM also have the same global accuracy because source zones that have their centroids in developed land are identical to those with developed land having the most area within them. For the BDAS method, the C-CAP data produced the best results, followed by the NLCD data, and then the CLEAR data. For the two IDM-based methods, C-CAP is the best ancillary data for the IDMC interpolator; CLEAR data has the lowest RMSE values for the IDMM interpolator, while NLCD data has the best results in terms of adj-RMSE, MAE, and adj-MAE. For the two GWR-based methods, either of which is usually the best among all methods, C-CAP results are the best among all LULC data sets for the GWRSC interpolator, while results based on the NLCD data are best for the GWRSP interpolator. For the three QR-based interpolators, C-CAP is the best ancillary data for the QRNC interpolator; however, the best data set varies with respect to the measurement scores for the other two QR-based estimators.

Table 1. Overall areal interpolation results derived from LULC data sets without enhancement

Interpolator	CLEAR				C-CAP				NLCD
	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE
BDAS	462.3	0.468	276.6	0.252	455.5	0.460	268.2	0.243	459.0	0.463	274.1	0.247
IDMC	461.0	0.466	275.8	0.251	450.9	0.448	269.0	0.242	469.5	0.473	281.5	0.255
IDMM	461.0	0.466	275.8	0.251	468.3	0.463	283.8	0.255	469.7	0.453	267.5	0.241
GWRSC	462.3	0.468	276.6	0.252	426.6	0.426	249.5	0.225	434.3	0.443	257.2	0.232
GWRSP	462.3	0.468	276.6	0.252	422.1	0.410	275.2	0.247	370.7	0.382	236.2	0.214
QRNC	462.3	0.468	276.6	0.252	440.3	0.435	262.6	0.237	450.1	0.456	273.7	0.247
QRSC	462.3	0.468	276.6	0.252	439.3	0.436	264.0	0.237	438.5	0.445	261.4	0.235
QRSP	462.3	0.468	276.6	0.252	485.5	0.436	304.7	0.272	450.1	0.449	284.9	0.258
Note: RMSE, root mean square error; Adj-RMSE, adjusted root mean square error; MAE, mean absolute error; Adj-MAE, adjusted absolute error; BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

Table 2 shows the global areal interpolation results derived from road buffer-enhanced LULC data sets. All the RMSE, adj-RMSE, MAE, and adj-MAE values presented in Table 2 are less than their corresponding values in Table 1. The reason that all GWR- and QR-based interpolators have the same results in this table is the same as that described for Table 1. The QRSP interpolator with C-CAP thematic ancillary data has the greatest absolute decrease of RMSE and MAE values, while the IDMM interpolator with C-CAP ancillary data has the greatest decrease in adj-RMSE values. IDMC with NLCD ancillary data has the greatest decrease of adj-MAE values. GWRSP using NLCD as ancillary data has the greatest decrease of RMSE and adj-RMSE values; QRSC with NLCD ancillary data has the least decrease of MAE values; and QRNC with C-CAP ancillary data has the least decrease of adj-MAE values.

Table 2. Overall areal interpolation results derived from LULC data sets enhanced by road buffer

Interpolator	CLEAR				C-CAP				NLCD
	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE
BDAS	399.7	0.398	235.2	0.212	393.4	0.391	235.1	0.213	394.0	0.395	226.9	0.207
IDMC	400.2	0.399	235.2	0.213	384.8	0.379	229.0	0.207	400.8	0.404	232.4	0.212
IDMM	400.2	0.399	235.2	0.213	396.5	0.386	248.5	0.225	382.3	0.381	223.8	0.203
GWRSC	399.7	0.398	235.2	0.212	375.4	0.364	227.1	0.204	373.7	0.370	218.2	0.198
GWRSP	399.7	0.398	235.2	0.212	382.8	0.366	237.1	0.211	358.8	0.355	209.2	0.188
QRNC	399.7	0.398	235.2	0.212	389.5	0.380	246.3	0.222	391.3	0.395	240.7	0.219
QRSC	399.7	0.398	235.2	0.212	390.3	0.379	244.5	0.219	389.6	0.394	255.9	0.216
QRSP	399.7	0.398	235.2	0.212	385.6	0.374	246.7	0.222	389.9	0.393	242.6	0.220
Note: RMSE, root mean square error; Adj-RMSE, adjusted root mean square error; MAE, mean absolute error; Adj-MAE, adjusted absolute error; BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

Finally, Table 3 shows the global areal interpolation results derived from parcel-enhanced LULC data sets. As expected, the overall error values in this table are the smallest compared to their corresponding values in the previous two tables. The GWRSC interpolator with the CLEAR ancillary data has the greatest decrease of RMSE, adj-RMSE, and MAE values, while the QRSP interpolator with the C-CAP ancillary data has the greatest decrease of adj-MAE values. IDMC with the NLCD ancillary data has the least decrease of RMSE, MAE, and adj-MAE values, while GWRSP with NLCD ancillary data has the least decrease in adj-RMSE values.

Table 3. Overall areal interpolation results derived from LULC data sets enhanced by parcel data

Interpolator	CLEAR				C-CAP				NLCD
	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE	RMSE	Adj-RMSE	MAE	Adj-MAE
BDAS	309.6	0.273	211.1	0.183	328.2	0.279	215.4	0.185	325.2	0.282	218.5	0.187
IDMC	307.0	0.258	214.5	0.182	338.0	0.274	222.1	0.188	390.1	0.309	271.6	0.227
IDMM	331.9	0.297	237.6	0.209	336.0	0.271	230.5	0.197	327.1	0.275	223.4	0.192
GWRSC	284.8	0.249	195.2	0.170	302.5	0.264	194.9	0.169	273.9	0.250	181.3	0.159
GWRSP	288.3	0.251	196.9	0.171	320.2	0.289	215.3	0.192	284.4	0.262	196.0	0.174
QRNC	315.6	0.282	223.4	0.196	320.0	0.276	206.9	0.180	317.8	0.273	213.1	0.184
QRSC	307.8	0.273	215.7	0.188	296.3	0.253	190.6	0.165	288.6	0.250	190.3	0.163
QRSP	309.1	0.275	216.2	0.189	326.8	0.277	210.2	0.181	334.4	0.284	218.9	0.187
Note: RMSE, root mean square error; Adj-RMSE, adjusted root mean square error; MAE, mean absolute error; Adj-MAE, adjusted absolute error; BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

The local evaluation of areal interpolation results was done by mapping the spatial distributions of absolute deviation errors for each target zone. Figure 2 compares the spatial distributions of the absolute errors of the areal interpolation results. The three maps in the first row of Figure 2 are the results derived from the three different remotely sensed data sets without any enhancement. The areal interpolation method with the lowest MAE value for each LULC data set in Table 1 was selected for mapping. The three maps in the second row are the results derived from road buffer-enhanced remotely sensed data. Again the areal interpolation methods for display were selected based on the same rule for values given in Table 2. Finally, the three maps in the last row compare the results derived from parcel-enhanced remotely sensed data. The selected criterion is also the smallest MAE value under each LULC data set. A four class, equal interval classification was used for each map. The visual complexity of each map in Figure 2 is very similar but the number of darker tones decreases from first row to last row.

Figure 2. The spatial distributions of absolute errors for areal interpolation results: (A) IDMC with CLEAR as ancillary data; (B) GWRSC with C-CAP as ancillary data; (C) GWRSP with NLCD as ancillary data; (D) BDAS with road buffer-enhanced CLEAR as ancillary data; (E) GWRSC with road buffer-enhanced C-CAP as ancillary data; (F) GWRSP with road buffer-enhanced NLCD as ancillary data; (G) GWRSC with parcel-enhanced CLEAR as ancillary data; (H) QRSC with parcel-enhanced C-CAP as ancillary data; (I) GWRSC with parcel-enhanced NLCD as ancillary data.

This visual result is verified by examining the statistical distributions of errors in each map class interval across all interpolators. Tables 4–6 present the statistical distributions of errors in each interval, and over- and under-predictions are further separated resulting in eight classes. Table 4 shows the statistical distributions of interpolation errors derived from remotely sensed data without any enhancement. Table 5 shows the statistical distributions of errors derived from road buffer-enhanced remotely sensed data. The general trend is that the observations move from tail intervals to the center intervals by applying the road buffer to CLEAR, C-CAP, and NLCD data. The results for the parcel-enhanced data sets remove even more observations from the classes at the tails. This helps to explain the greater reduction in RMSE values with enhancement than for the MAE values in Tables 1–3.

Table 4. The number of observations in each break interval of error maps derived from three different remotely sensed data sets

		<–300	−300 to –201	−200 to –101	−100 to 0	1–100	101–200	201–300	>300
CLEAR	BDAS	48	29	39	60	45	26	28	48
	IDMC	50	28	37	62	46	25	25	50
	IDMM	50	28	37	62	46	25	25	50
	GWRSC	48	29	39	60	45	26	28	48
	GWRSP	48	29	39	60	45	26	28	48
	QRNC	48	29	39	60	45	26	28	48
	QRSC	48	29	39	60	45	26	28	48
	QRSP	48	29	39	60	45	26	28	48
C-CAP	BDAS	44	29	42	56	46	41	22	43
	IDMC	52	29	29	67	43	39	15	49
	IDMM	59	22	35	66	35	38	17	51
	GWRSC	45	31	32	64	49	36	25	41
	GWRSP	56	25	35	54	39	39	23	52
	QRNC	48	29	36	61	51	28	20	50
	QRSC	51	27	37	60	44	38	18	48
	QRSP	59	25	36	54	37	36	23	53
NLCD	BDAS	49	27	48	46	48	29	27	49
	IDMC	48	30	49	51	38	34	21	52
	IDMM	51	24	45	55	42	37	19	50
	GWRSC	50	19	52	52	44	38	22	46
	GWRSP	44	22	45	51	54	42	23	42
	QRNC	52	27	37	59	42	39	13	54
	QRSC	49	23	49	48	49	37	21	47
	QRSP	56	25	33	52	59	26	17	55
Note: BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

Table 5. The number of observations in each break interval of error maps derived from three different remotely sensed data sets enhanced by road buffer

		<–300	−300 to –201	−200 to –101	−100 to 0	1–100	101–200	201–300	>300
CLEAR	BDAS	36	22	50	61	52	35	30	37
	IDMC	36	24	48	61	54	34	28	38
	IDMM	36	24	48	61	54	34	28	38
	GWRSC	36	22	50	61	52	35	30	37
	GWRSP	36	22	50	61	52	35	30	37
	QRNC	36	22	50	61	52	35	30	37
	QRSC	36	22	50	61	52	35	30	37
	QRSP	36	22	50	61	52	35	30	37
C-CAP	BDAS	38	26	49	60	41	41	30	38
	IDMC	40	22	48	59	44	44	30	36
	IDMM	45	23	52	52	36	48	23	44
	GWRSC	40	23	47	59	39	49	29	37
	GWRSP	41	29	47	50	40	43	31	42
	QRNC	46	27	40	59	34	44	31	42
	QRSC	45	28	41	55	36	45	31	42
	QRSP	45	29	43	55	36	40	33	42
NLCD	BDAS	35	27	38	68	55	31	33	36
	IDMC	37	30	41	66	47	35	31	36
	IDMM	39	18	50	66	41	45	23	41
	GWRSC	37	21	51	58	57	40	24	35
	GWRSP	34	23	46	68	56	39	25	32
	QRNC	40	32	40	54	48	41	24	44
	QRSC	39	32	44	52	44	45	28	39
	QRSP	44	29	42	53	42	43	27	43
Note: BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

Table 6. The number of observations in each break interval of error maps derived from three different remotely sensed data sets enhanced by parcel.

		<–300	−300 to –201	−200 to –101	−100 to 0	1–100	101–200	201–300	>300
CLEAR	BDAS	40	27	33	59	62	36	30	36
	IDMC	39	28	41	53	68	30	23	41
	IDMM	51	26	28	66	52	30	26	44
	GWRSC	32	22	48	69	62	34	24	32
	GWRSP	33	22	48	68	61	34	25	32
	QRNC	42	25	44	73	38	38	19	44
	QRSC	37	27	47	67	46	37	22	40
	QRSP	38	26	46	68	46	37	22	40
C-CAP	BDAS	38	24	41	59	49	54	25	33
	IDMC	37	32	38	61	54	37	31	33
	IDMM	37	29	44	57	48	42	28	38
	GWRSC	33	25	39	62	61	49	26	28
	GWRSP	41	27	35	63	54	39	34	30
	QRNC	36	27	36	66	68	32	23	35
	QRSC	38	14	42	64	68	44	26	27
	QRSP	33	29	39	60	67	38	24	33
NLCD	BDAS	40	29	25	57	69	37	30	36
	IDMC	56	22	42	54	42	38	20	49
	IDMM	38	27	41	59	61	31	32	34
	GWRSC	27	27	35	76	65	48	19	26
	GWRSP	31	27	49	63	60	36	26	31
	QRNC	32	28	44	70	50	45	18	36
	QRSC	34	22	40	65	67	42	22	31
	QRSP	34	30	39	73	55	36	19	37
Note: BDAS, binary dasymetric; IDMC, intelligent dasymetric mapping centroid; IDMM, intelligent dasymetric mapping most; GWRSC, geographically weighted regression with scaled coefficients; GWRSP, geographically weighted regression with scaled population; QRNC, quantile regression with non-negative coefficients; QRSC, quantile regression with scaled coefficients; QRSP, quantile regression with scaled population.

As for the results derived from these three data sets, most numbers under the two tail break intervals in Table 4 are larger than the corresponding ones in Table 5, while most numbers under the two center break intervals in Table 4 are less than the corresponding ones in Table 5. However, this trend does not occur when applying the road buffer to the C-CAP data. Although the observations in the two tail intervals decrease, the observations within the two center intervals also decrease when applying the road buffer to the C-CAP data. In terms of the statistical distributions of interpolated population errors, the road buffer enhancement has a greater effect on CLEAR and NLCD data than on C-CAP data.

6. Conclusions

Enormous effort and cost is spent compiling publicly available LULC data set from remotely sensed imagery. When deciding on which categories to include in the classification, it is important to anticipate the uses of these data sets in further analysis. Areal interpolation is a widely used method for rectifying disparate data sets in GIS. This article evaluates the effect of three publicly available land cover data sets on intelligent areal interpolation results and, at the same time, evaluates the effects of two enhancements on areal interpolation results by applying them to remotely sensed land cover data. Results show that national and regional data sets can be of greater use in areal interpolation if more attention is given to the division of developed urban land into finer categories versus a single broad category.

Based solely on the original classified data sets, NLCD and C-CAP produced better results than CLEAR because CLEAR did not provide any sub-categorization of developed urban land. However, the advantage of C-CAP and NLCD data over CLEAR data becomes less obvious if these data sets were enhanced using additional information. The enhancement of CLEAR data with a parcel layer improved CLEAR results, making it better than C-CAP for six of the eight methods and better than NLCD for four of the eight methods with respect to RMSE. CLEAR is only better than the national data set for three methods based on MAE.

Overall, using parcels appears to be a better spatial enhancement than a road buffer. One reason is that road buffer enhancement constrains only the developed land spatially, while the parcel enhancement not only constrains the developed land spatially, but also expands the developed land into more categories. The developed land category in CLEAR, for example, was further divided into single-family developed, two-family developed, and multiple-family developed. Another reason for parcel data being a better enhancer of remotely sensed land cover data is that the developed land within the road buffer area might not contain only residential land along the road, but also the road itself.

Ancillary data are critical to improving areal interpolation results and are as important as the actual interpolation method used. In this study, overall weaker methods, such as BDAS, the two IDM methods, and QRSP outperformed the overall stronger methods when the original remotely sensed data were enhanced by road buffers or parcel data layers. For areal interpolation, having ancillary data that are spatially correlated variables is very important. More effort is needed to improve the spatial detail among developed urban categories. Parcel data provide a better enhancement than road buffers because the parcels remove the areas of the roads themselves. Land cover data sets based on the broadest Anderson Level I classification such as the CLEAR data set in this study are insufficient for use as ancillary data to estimate population distributions. Enhancement before using imagery data in an areal interpolation analysis, which can be implemented easily in a GIS, is a solution to this problem.