Transactions in GIS, 2010, 14(3): 351–377
Research Article
ASPHAA: A GIS-Based Algorithm to Calculate Cell Area on a LatitudeLongitude (Geographic) Regular Grid Monia Santini
Andrea Taramelli
Department of Forest Environment and Resources University of Tuscia Viterbo, Italy
ISPRA – High Institute for Environmental Research Roma, Italy
Alessandro Sorichetta Dipartimento di Scienze della Terra “Ardito Desio” Università degli Studi di Milano Milano, Italy
Abstract One characteristic of a Geographic Information System (GIS) is that it addresses the necessity to handle a large amount of data at multiple scales. Lands span over an area greater than 15 million km2 all over the globe and information types are highly variable. In addition, multi-scale analyses involve both spatial and temporal integration of datasets deriving from different sources. The currently worldwide used system of latitude and longitude coordinates could avoid limitations in data use due to biases and approximations. In this article a fast and reliable algorithm implemented in Arc Macro Language (AML) is presented to provide an automatic computation of the surface area of the cells in a regularly spaced longitude-latitude (geographic) grid at different resolutions. The approach is based on the well-known approximation of the spheroidal Earth’s surface to the authalic (i.e. equal-area) sphere. After verifying the algorithm’s strength by comparison with a numerical solution for the reference spheroidal model, specific case studies are introduced to evaluate the differences when switching from geographic to projected coordinate systems. This is done at different resolutions and using different formulations to calculate cell areas. Even if the percentage differences are low, they become relevant when reported in absolute terms (hectares). tgis_1200
351..378
Address for correspondence: Andrea Taramelli, ISPRA Institute for Environmental Protection and Research, via di Casalotti, 300, 0166, Roma, Italy. E-mail:
[email protected] © 2010 Blackwell Publishing Ltd doi: 10.1111/j.1467-9671.2010.01200.x
352
M Santini, A Taramelli and A Sorichetta
1 Introduction Global lands span over an area greater than 15 million km2 and information types are highly variable. A variety of protocols and conventions have been established in order to deal with comprehensive aspects of the environment (Mitchell 2003). An increasing number of anthropogenic and biophysical phenomena occurring at multiple spatial and temporal scales, such as land cover/use transformations or climate changes, as well as their assessment for atmospheric studies or hazard-vulnerability-risk analyses, are treated as an integrated system since years (United Nations 1951, Gringorten 1972, Henderson-Sellers et al. 1986, Kimerling et al. 1999, McCallum et al. 2006). To this regard the Geographic Information System (GIS) is very helpful as it both addresses the necessity to handle a large amount of data at multiple scales and allows the spatio-temporal integration of datasets deriving from different sources. Modeling, in particular, benefits from GIS for either assessment studies or future scenarios about environmental problems, useful in planning prevention, mitigation and adaptation strategies. As an example, the application of numerical prediction models to strong weather phenomena such as dust and sand storms is considered of key importance in developing control/mitigation measures (Engelstaedter et al. 2006, Washington et al. 2006, Bach et al. 2007) in regions far from the source areas. These and many others models, for example the General or Regional Circulation Models (Pasqui et al. 2004, 2005, FoxRabinovitz et al. 2006) or the Dynamic Global Vegetation Models (DGVMs; e.g. Sitch et al. 2003, Crevoisier et al. 2007) run on latitude-longitude (geographic) grids over which the areal extent of model inputs and parameters needs to be carefully estimated to provide more reliable results. In other cases, principally for hazard/risk mapping and predictions at scales from global to regional, it is reasonable to work on a geographic basis given the size of involved areas (Dilley et al. 2005). An accurate computation of the cell area on a geographic grid-based map is necessary to appropriately assess the extent of critical conditions up to detecting hazard/risk hot-spot distribution (Dilley et al. 2005, Arnold et al. 2006). All these studies have a common baseline: interpolating the raw observations, the post-processed data and the model inputs and outputs into a shared grid-structure (Snyder 1926, Meshcheryakov 1962, Kimerling et al. 1999, Sanyal and Lu 2004, Yuan 2005). Using geographic coordinate systems, many area-based socio-economic and environmental indicators represented in grid format (e.g. population density, biomass consumption, crop yield, pollutant dispersal) can be evaluated for global analyses with higher accuracy. The issue of transforming the coordinate system of a grid, whose reference change and planimetric projections are critical steps as source of errors and approximations, suggests that geographic coordinates should be preferred to calculate area, even in place of equal-area projections (Yetman et al. 2000, Small and Cohen 2004, Hay et al. 2006, McCallum et al. 2006, Elvidge et al. 2007). To this effort a global geographic grid is useful for partitioning the Earth into a geometrically regular set of cells (Tobler and Chen 1986, Bonham-Carter 1994, Burrough 1997, 2001), after selecting the most suitable Earth’s shape representation. Although the geoid (i.e. the irregular surface perpendicular at each point to the force of gravity) is the solid better approximating the elevation contours of the globe, most GIS © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
353
databases, used for cartographic-planimetric purposes, rely on ellipsoidal (regular) models of the Earth. In particular, such ellipsoids are simplified to oblate spheroids, whose two horizontal axes have the same length and are longer than the third vertical axis (Thomas 1952). The reference surface of the spheroid, known as the “Datum”, varies to best fit the geoid at different Earth’s locations and it is characterized by a fixed initial point (origin), the distance between the geoid and the Ellipsoid at such a point and the initial azimuth on the surface. A further simplification allows approximating the Earth’s spheroid to the related authalic sphere, being a perfect sphere having the same surface area as the Earth’s spheroid (Kimerling et al. 1999). Each portion of the spheroid surface extending between two meridians and two parallels maintains the surface area in its corresponding “projected image” on the authalic sphere. The possibility to choose between spherical or spheroidal models introduces a key issue. Latitude coordinates provided by Global Position Systems (GPS) and stored into geographic databases are usually defined “geodetic” (or “geographic”), meaning that the latitude angle from the equator is measured from the intersection point between the equatorial plane and the line crossing normally the spheroid surface at the considered point. The so-called “geocentric” latitude is instead measured from the center of the Earth and is more suited for sphere-based approaches. For each geodetic latitude on the spheroid surface the authalic latitude on the simplified spherical surface can be computed. Hereafter the spheroid model will be indicated as SPHD while the related authalic sphere as ASPH. The ASPH can be partitioned into a geographic regular grid (Figure 1a), partially matching the criteria reported in Kimerling et al. (1999), whose cells cover the entire globe without overlapping, and with the same topology and shape. Each portion of the Earth’s surface can be represented by cells having the same angular dimensions along the NS and EW directions (i.e. spherical quadrangles or triangles where meridians converge to the poles). Equally-spaced latitude degree intervals of the regular geographic grid represent the rows of the grid, while the equally-spaced longitude intervals represent the columns of the grid. In this representation the width of the longitude is proportional to the cosine of the latitude, while the length of a given angle of latitude remains constant (equal to the authalic sphere radius multiplied by the cell dimensions in radians; Figure 1b). This represents a well-known issue to solve when using a latitude-longitude geographic grid, as grid cells decrease in area towards the poles. Although it was widely demonstrated that such a geographic grid provides nearly error-free conditions for area estimation (Snyder 1926, Perillo et al. 1999), the computational procedure to define the area of the grid cells is crucial, as the units of measure of cell size and cell surface area are not the same. In this article a fast and reliable algorithm implemented in Arc Macro Language (AML; ©ESRI, Redlands, California) syntax is presented to provide the automatic computation of the cell area in a regularly spaced geographic grid at multiple spatial scales. Such an approach is based on the above mentioned approximation of spheroid Earth’s surface to the authalic sphere. In the next sections the algorithm is described, tested and applied to key case studies in order to highlight how areas estimated from grids in Projected Coordinate Systems (PCSs) suffer from biases in comparison with the ones computed from grids in Geographic Coordinate Systems (GCSs), and that such deviations differ according to latitude, © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
354
M Santini, A Taramelli and A Sorichetta a
b
Figure 1 (a) Example of the partitioning of a sphere surface into a grid with a cell size dimension equal to 10° along both the X (longitude) and the Y (latitude) dimensions. Where meridians converge to the poles, the cells are spherical triangles; and (b) on a geographic grid meridians converge toward the poles, whereas the parallel interdistance remains constant (modified from http://members.shaw.ca/quadibloc/maps/ mcy0105.htm)
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
355
Figure 2 Example of a generic cell of the grid comprised between l1 and l2 degrees of longitude and between j1 and j2 degrees of (authalic) latitude. l and j are generic coordinates of a point belonging to the cell, whereas R is the (authalic) Earth’s radius (modified from http://badc.nerc.ac.uk/help/coordinates/cell-surf-area.html)
cell resolution and projection scheme. Finally the differences between ASPHAA and a different cell area computation approach are evaluated.
2 Algorithm Description The surface area S of a generic cell of the aforementioned grid, comprised between l1 and l2 degrees of longitude and between j1 and j2 degrees of (authalic) latitude (i.e. the green cell in Figure 2) can be calculated as (Tobler and Chen 1986): ϕ2
S=
ϕ
ϕ
λ2 2 2 ⎛ λ2 ⎞ 2 2 ( ) R R cos ϕ d λ d ϕ = R cos ϕ d ϕ d λ = R λ − λ 2 1 ∫ cos ϕ dϕ ⎟ ∫ ⎜⎝ λ∫ ∫ ∫ ⎠ ϕ1 ϕ1 λ1 ϕ1 1
(1)
giving:
S = R 2 ( λ2 − λ1 ) (sin ϕ 2 − sin ϕ1 ) where
0 ⱕ l1 ⱕ l ⱕ l2 ⱕ 2p -p/2 ⱕ j1 ⱕ j ⱕ j2 ⱕ p/2 R = the (authalic) radius of the sphere j1, j2 and l1, l2 are expressed in radians (p radians = 180 degrees). © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
(2)
356
M Santini, A Taramelli and A Sorichetta a
b
Figure 3 (a) This figure represents an equal angle grid having the same X and Y spacing (CS) in angular dimensions. For each row, the coordinates of the bottom edge are computed from the coordinates of the bottom of the grid (YLL) plus the value (NROWS - NROW) times the cellsize (CS). The coordinates of the top edge are computed by YLL + (NROWS - NROW + 1) * CS; and (b) the $$rowmap is a built-in grid in ArcInfo that can be used to identify row locations of a grid cell. The $$rowmap grid assigns 0 to all the cells in the first (upper) row, 1 to the cells in the second row, and so forth
Based on the above assumption it is possible to verify that all the cells belonging to the same row of the grid (i.e. comprised in the same interval of latitude) have the same surface area: the three factors of Equation (2) (R2, (sin j2 - sin j1) and (l2 - l1), respectively) are constant. Thus, the grid area computation can be simplified as it only requires the evaluation of the cell surface associated to a single column of the grid, being that such an estimate is valid for all the other columns. Within one grid column, the surface area for each cell can be calculated starting from its upper and lower (authalic) latitudes. The whole procedure is better described in the following through an explicative example. The example considers a grid with five rows and seven columns (Figure 3a), the variable NROW indicates the row numbering on the grid (from top to bottom), NROWS is the total number of rows (i.e. 5), CS is the cell size (in radians), YLL is the (authalic) latitude coordinate of the lower left corner of the grid (in radians), and R is the (authalic) Earth’s radius in m. When using GIS databases, before applying Equation (2) it is often © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
357
necessary to convert the geodetic coordinates stored by the grid reference system into authalic coordinates. The cell area for the lowest row of the grid (the dark grey one in Figure 3a), considering Equation (2), is:
R ⋅ R ⋅ CS ⋅ [sin ( YLL + 1⋅ (CS)) − sin ( YLL )]
(3)
while the cell area for the next row (the light grey one) of the grid, still from Equation (2), is:
R ⋅ R ⋅ CS ⋅ [sin ( YLL + 2 ⋅ (CS)) − sin ( YLL + 1⋅ (CS))]
(4)
Such a method is coded in an AML script (Appendix A1) in order to derive, from a generic latitude-longitude grid supplied in GRID format (ESRI 1990), a grid storing in each cell the area for that same cell in square linear units (m2). The procedure can be divided into six steps as follows. 1. The pi-Greco (p) value is fixed to 3.14159265358979, the CS and the Geodetic Lower Left (GLL) corner coordinates are derived from the grid metadata and converted into radians and the total number of rows (NROWS) is stored 2. The native geodetic Datum related to the dataset is selected among 24 coordinate systems. For each one, the semi-major axis (SMA) and the flattening (FL) of the spheroid are considered 3. From FL the eccentricity e is computed as
e = ( 2 ⋅ FL ) − (FL2 )
(5)
4. A further grid named ROW is created, spatially matching the input grid and keeping, for each cell, the value of the row occupied by that same cell (counting the rows from top to bottom). This is achieved by computing
ROW = $$rowmap + 1
(6)
where $$rowmap is an ArcInfo built-in grid that can be used to identify the row numbering, but as it starts from 0 (Figure 3b), one unit is added to the whole grid in order to get the right numbering starting from 1. 5. Recalling that cells belonging to the same row also have the same surface area, the authalic latitudes for the lower and upper bound of each row (j1 and j2, respectively) are computed from (Snyder 1926):
ϕ1 = arcsin ( q1 q )
(7)
ϕ 2 = arcsin ( q 2 q )
(8)
sin ϕ g 1 1 ⎤ ⎡ (1 − e ⋅ sin ϕ g1 ) ⎤ ⎡ q1 = (1 − e2 ) ⋅ ⎢ − ⎥ ⋅ ln ⎢ ⎥ 2 ⎣ 1 − e ⋅ sin 2ϕ g 1 2e ⎦ ⎣ (1 + e ⋅ sin ϕ g1 ) ⎦
(9)
where
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
358
M Santini, A Taramelli and A Sorichetta sin ϕ g 2 1 ⎤ ⎡ (1 − e ⋅ sin ϕ g 2 ) ⎤ ⎡ q 2 = (1 − e2 ) ⋅ ⎢ − ⎥ ⋅ ln ⎢ 2 1 − e ⋅ 2 2 e ⎦ ⎣ (1 + e ⋅ sin ϕ g 2 ) ⎥⎦ sin ϕ g2 ⎣
()
( (
⎡ sin π ⎤ ⎡ 1 − e ⋅ sin 1⎥ ⎢ 2 2 ⎢ q = (1 − e ) ⋅ ⎢ − ⎥ ⋅ ln ⎢ 2 ⎢ 1 − e ⋅ sin π 2e ⎥ ⎢ 1 + e ⋅ sin ⎣ ⎦ ⎣
( π )) ( π ))
⎤ 2 ⎥ ⎥ ⎥ 2 ⎦
(10)
(11)
where jg1 and jg2 are lower and upper geodetic latitudes, respectively. In the case of the lowest row of the grid, jg1 ≡ GLL. 6. From SMA and q the authalic radius R is computed as (Snyder 1926):
R = SMA 2 ⋅
q 2
(12)
Now equation (2) can be applied. According to the ArcInfo specifications, the upper limit of the grid size to run this algorithm is assumed to be 4,000,000 of rows and 4,000,000 of columns. This would mean a grid working at ca. 0.324″ lon ¥ 0.162″ lat cell resolution (ca. few meters) for the whole globe, which would be computationally a very hard challenge for global applications. However if the angular spacing (resolution) of grid cells is so large that areas computed in m2 overcome the 32-bit limit of representation allowed for grids in ArcInfo, the computation can be also made in km2 by converting the SMA values from m to km.
3 Algorithm Testing Even if the reliability of the ASPH in approximating the SPHD surface areas is widely assessed (e.g. Tobler et al. 1995, Qihe et al. 1999), the amount of biases in estimating the cell surface area as a function of the grid resolution and latitude did not receive great attention. Before testing the ASPHAA, a comparison was made between its results and a numerical solution of the grid-cell area taking as reference the source SPHD. A Fortran 90 routine for the computation of the cell area on a regular (geographic) grid partitioning a SPHD into cells is implemented. The aim of this procedure was: (1) defining the performances of the algorithm for different cell resolutions and latitudes in order to assess its accuracy when used for different study areas and scales; and (2) as the algorithm is working on ArcInfo format, which is able to store grid values with a 32-bit (float) precision, the biases in case of large cells (i.e. billion of sq. meters) have to be taken into account for uncertainty evaluation. As it was already highlighted for the ASPH, also for the SPHD all the columns of the regular grid have the same surface area. There is also a perfect symmetry between upper (north) and lower (south) semi-SPHDs. For this reason and in order to meet computational effort requirements, only a single longitude column lying from 0° to 90°N latitude (semi-lune) was considered (Figure 4a). Whatever is the chosen angular grid resolution, such a numerical procedure divides each SPHD’s cell into further equal angle cells. Aiming to a precision of 1 m2 for area © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm a
359
b
Figure 4 (a) In grey an example of semi-lune on a sphere; and (b) example of spherical triangle
computation, it was first acknowledged that a 1″ ¥ 1″ cell represents a spherical trapezium from 0° to 89°59′59″N latitude and a spherical triangle from 89°59′59″ to 90°00′00″ N latitude (the northern side is one point, i.e. the pole; Figure 4b). At similar scales (1″ ¥ 1″ is well representing a 30 m ¥ 30 m square at the equator) it is possible considering a purely 2D cell for which the biases from the real Earth’s curvature can be neglected. Thus the chosen numerical solution uses as minimum partitioning unit 1″ ¥ 1″ cells (arc-cells hereafter). The routine reads the grid resolution, the lower left corner geodetic latitude and longitude and, for each cell along the single grid column, it first calculates the area for each arc-cell it contains and then for the cell itself. The algorithm’s pseudo-code is reported in Appendix A2. Seven grid resolutions were selected: 0.1, 0.25, 0.5, 1, 2.5, 5 and 10°, whose cells contain, respectively, 360 ¥ 360, 900 ¥ 900, 1,800 ¥ 1,800, 3,600 ¥ 3,600, 9,000 ¥ 9,000, 18,000 ¥ 18,000 and 36,000 ¥ 36,000 arc-cells. For each resolution the total surfaces both of the “numeric” SPHD and the ASPH were calculated by multiplying the semi-lune area by the appropriate factor (e.g. 360 ¥ 2 in the case of 1° of resolution, 720 ¥ 2 for 0.5°, etc.). Both calculations show a good agreement (deviations less than 10-5–10-6%) from the SPHD area analytically calculated as:
⎧ ⎡ ⎪ ⎢ (SMI 2 ) A = 2π ⋅ ⎨SMA 2 + ⎢ SMI ⎪ ⎢ sin arccos ⎩ ⎣ SMA
(
( ))
(
( )) ( )
⎤ ⎫ ⎧ ⎡1 + sin arccos SMI ⎤ ⎫ ⎥ ⎪ ⎪ ⎣⎢ SMA ⎦⎥ ⎪ ⎬ ⎥ ⎬ ⋅ ln ⎨ SMI ⎪ ⎥⎪ ⎪ ⎦⎭ ⎩ ⎭ SMA
(13)
where SMA and SMI are semi-major and the semi-minor axes, respectively. Then, taking the SPHD model as reference, for each resolution the mean Absolute Percent Deviation (APD) and Root Mean Square Error (RMSE) of ASPH were computed. Graphs in Figure 5 report exponential curves that fit the relationships (blue points) of the grid resolution with APD (R2 = 0.989) and RMSE (R2 = 0.988). In the APD case it is © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
360
M Santini, A Taramelli and A Sorichetta
Figure 5 In blue are reported the APD (left) and RMSE (right) of ASPHAA algorithm from SPHD-based numerical solution of the cell area for different angular cell resolutions. The red point represents the same evaluations made in the work of Tobler et al. (1995). The fitted line refers only to the blue points
noteworthy that when the angular cell spacing increases the biases among the two – ASPH-based and SPHD-based – algorithms decrease. For the RMSE, as expected, the behavior is inverted: the minor percent deviation at coarser resolutions corresponds to a much larger reference area. Absolute ASPHAA deviations appear raised when latitude increases (not showed here), this being more evident at finer resolutions, and suggesting that resolution- and latitude-related correction factors could be introduced to improve ASPHAA evaluations. It is easy to extrapolate that in the case of a resolution finer that 0.1° the above trends are respected. A concrete example is reported in the Global Demography Project made for the year 1995 by Tobler et al. (1995). Starting from the grid storing the people living in each cell at 5′ ¥ 5′ (ca. 0.083° ¥ 0.083°) of resolution, the population density grid was calculated by dividing the abovementioned population count grid by the cell area grid, and considering an authalic spherical surface with a radius of 6,731,007.178 m. This result was compared with the one derived dividing the same population count grid by the cell-area grid calculated with ASPHAA. The maximum difference consisted of one inhabitant and occurred in 0.23% of cells. In Figure 5 the APD and RMSE for the sphereand spheroid-based algorithms are reported (red points). When included in the fitting curves previously computed for APD and RMSE, their R2 decreases to 0.896 and 0.927, respectively. It is interesting to note that authalic latitudes were not used in the case of sphere-based computations of Tobler et al. (1995).
4 ASPHAA Applications In this section several application examples are shown in order to assess when and how the GIS-based ASPHAA could be useful within different research topics, overcoming the limits of both computing grid cell area by other procedures and using databases in projected coordinate systems. In particular, algorithm evaluation is given in two directions: first, there is an assessment of its performance at different spatial scales (global, continental and regional) and second, the influence of area calculation procedures on the © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
361
distribution assessment of biophysical and socio-economic factors, considered both separately and coupled, and commonly used for risk analyses, is evaluated. The aim is to highlight the impact on area calculation due to the data source type and/or aggregation.
4.1 Global Applications 4.1.1 Crop production In the context of global change studies, modifications in land use pressures were recently recognized to have a key role in greenhouse gas emission increase or reduction (McGuire et al. 2001, Albani et al. 2006). One of the main research topics is the tentative coupling between biophysical aspects of land use changes and their socio-economic drivers. This has been done, for example, by Ronneberger et al. (2009), who tried to integrate the global agricultural Kleines Land Use Model (KLUM; Ronneberger et al. 2005) and the Global Trade Analysis Project (GTAP) model (Hertel 1997). A useful dataset for similar evaluations was built by Monfreda et al. (2008), which combined the FAOSTAT (FAO 2006a) and Agro-MAPS (FAO 2006b) databases with a number of national censuses and surveys to create an extensive global dataset on crops for the year 2000. Each cell of their grid-based dataset contains the fractions of harvested area (hectares) and the yield (metric tons per harvested hectare) per year for 175 individual crops at a 5′ ¥ 5′ spatial resolution (http://www.geog.mcgill.ca/~nramankutty/Datasets/Datasets.html). In some cells a crop is harvested multiple times in a single year so that its surface is counted more than once and its cell fractions can exceed one. While the crop harvested area could be useful as an initial condition for most of land use change models, the crop yield attribute is probably most helpful for socio-economic evaluation driving land use demands through market dynamics. Since such a spatial dataset was derived by merging census statistics from different sources and scales and implementing spatial analysis procedures to interpolate data and to fill voids (Monfreda et al. 2008), a direct comparison between aggregate (administrative level) and spatialized statistics is not easy to accomplish. The four globally most diffused crops, i.e. wheat, rice, maize and soybean (Monfreda et al. 2008), were chosen to compare ASPHAA-based areas with those resulting from projecting the grids with an equal-area global projection (Goode Homolosine, see Appendix A3 for projection parameters), maintaining the source Datum (WGS 1984) and choosing a 10,000 m mesh size resolution, which is comparable with 5′ cell angular dimensions. On average, the projected grid-based area shows a bias (overestimation) compared to the ASPHAA reference area of about 0.2% (ca. 1,129,000 ha). The yield is affected on the order of more than 3,800,000 tons being overestimated each year. Concentrating on the two most diffuse crops, i.e. wheat and rice, is it possible to observe that the former spans a much huger latitude interval (Figure 6a) than the latter (Figure 6b); this is the reason why it appears better when balancing the overestimations and underestimations by ASPHAA, resulting in lower percent differences with census data. Then, a test of the differences among areas (in decimal million of hectares) resulting from census databases and the ASPHAA algorithm is carried out using the aggregation of the 175 crops into 11 macro-categories as in Table 3 from Monfreda et al. (2008). Results (Figure 7) show that, on average, the ASPHAA algorithm underestimates the area for 1%, and for each crop group the percentage increase with the total surface extent slightly following a logarithmic law (R2 = 0.68). In evaluating such differences one has to © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
362
M Santini, A Taramelli and A Sorichetta
a
b
Figure 6 Map of harvested areas (m2) per year related to (a) wheat and (b) rice
consider that biases can be due to the transformation of administrative unit boundary maps from vector to grid and to the interpolation procedures that propagate the error of merging crops into macro-categories.
4.2 Risk Analysis A further interesting research area, which integrates biophysical and human factors, is the analysis of risk related to natural hazards. After the Indian Ocean Tsunami in 2004 and the South-east Asia earthquakes and Central America hurricanes in 2005 (Arnold et al. 2006), the Global Natural Disaster Risk Hotspots (GNDRH) project was launched by the World Bank and Columbia University Center for International Earth Science © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
363
Figure 7 Percent differences between areas (Mha) reported in Monfreda et al. (2008) and computed by ASPHAA after aggregating 175 crops into 11 macro-categories
Information Network (CIESIN). Here a noticeable effort was made in order to produce a comprehensive geodataset attempting to assess multi-hazard risk analysis at a globalscale and to identify the hot-spots, which are regions where the risk of natural disasters is particularly high or regions would impact other adjacent areas. Data are presented in the form of gridded global maps with a 2.5′ ¥ 2.5′ spatial resolution (http:// www.ldeo.columbia.edu/chrr/research/hotspots/coredata.html). The risk level was evaluated by combining the natural hazard probability and the historical vulnerability for the considered area with the exposure of elements at risk (lives, properties, etc.). In particular the global risk assessment focused on two disaster consequences, mortality and economic losses, for six major natural hazards: earthquakes, volcanoes, landslides, floods, droughts and cyclones. In order to provide a relative representation of risk, absolute numbers (hazard frequency, mortality and economic losses) were converted into index values from 1 to 10 corresponding to deciles of their distribution. Finally a multi-hazard risk index was created by summing decile category values among all six natural hazards. This multi-hazard index reflects the number of hazards assumed relatively significant in a particular grid cell (Dilley et al. 2005). Considering as attributes both the frequency of the six natural hazards and their risk in terms of mortality and economic losses, areas belonging to each attribute decile and for each type of hazard were calculated both with ASPHAA and after projecting geographic grids into an equal-area world projection (i.e. Sinusoidal; see Appendix A3 for projection parameters), maintaining the source Datum (i.e. WGS 1984) and choosing a 5,000 m mesh size resolution, comparable with 2.5′ cell angular dimensions. Percent deviations of the projected- from the ASPHAA-based areas are reported in the graphs of Figures 8a–c showing significant differences mainly in the evaluation of hazard consequences (economic losses and mortality), lower for floods, droughts and earthquakes and higher for cyclones, landslides and volcanoes. When considering multi-hazards indices, the difference in percent seems to increase with the increase of the multi-hazard index itself for frequency and mortality, following an exponential law and having correlation coefficients equal to 0.78 © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
364
M Santini, A Taramelli and A Sorichetta a
b
c
Figure 8 Percent differences between areas computed using PCS and GCS for the different risk deciles regarding hazard: (a) frequency; (b) mortality; and (c) economic losses
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
365
and 0.74, respectively. For economic losses, it is more difficult to detect a trend since numerous risk classes disappear when switching from geographic to a projected coordinate system.
4.3 Continental Applications 4.3.1 Dynamical global vegetation models The first case at a continental scale considers a widely used Dynamical Global Vegetation Model (DGVM), the Lund-Postdam-Jena (LPJ) scheme (Sitch et al. 2003). LPJ is a process-based model simulating the water and carbon exchanges between the biosphere and atmosphere by means of a given set of parameters and input variables. For each grid cell, vegetation is described in terms of the fractional coverage of 10 different Plant Functional Types (PFTs) that are able to compete for space and resources. The PFTs may in principle co-exist at any location, depending on plant competition and a set of environmental constraints. Their relative proportion is determined by the competition among vegetation types with specific ecological strategies for dealing with temperature, water and light. Among numerous land classification systems, the MODerate resolution Imaging Spectroradiometer (MODIS) sensor data elaborations supplied by NASA (http:// modis.gsfc.nasa.gov/data/dataprod/index.php), offer the land cover classification useful to represent LPJ’s PFTs. As a proxy of a typical PFT partitioning reproduced through LPJ simulations, five maps of fractional cover (referring to Crop, C; Grassland and Wetlands, GW; Evergreen Broadleaved Forests, EBF; Evergreen Needleleaved Forests, ENF; and Deciduous Broadleaved Forests, DBF) were here used for Europe in a grid format and at a 0.25° ¥ 0.25° spatial resolution. First for each pixel and then for the whole map, the total area (in m2) covered by the five vegetation classes was computed, either using ASPHAA or the approach implemented in LPJ, which calculates the area for a pixel on a lat/lon grid as (http://www.pikpotsdam.de/research/cooperations/lpjweb):
S = (111, 000 ⋅ CS) ⋅ (111, 000 ⋅ CS) ⋅ cos ( LATc )
(14)
where LATc is the geodetic latitude of the cell center. Table 1 shows the differences in absolute area and percent between the two algorithms. It is worth noting the cumulated underestimation of the LPJ approach for all the
Table 1 Absolute and percent differences between total surfaces computed by the procedure implemented in LPJ and the ASPHAA algorithm PFT
ASPHAA (m2)
LPJ’s algorithm (m2)
Diff (km2)
Diff (%)
C DBF EBF ENF GW TOT
4,819,400,523,776 1,329,023,811,584 65,319,903,232 2,870,795,304,960 1,428,607,336,448 10,513,146,880,000
4,788,854,980,608 1,319,833,960,448 64,964,202,496 2,846,025,580,544 1,420,176,785,408 10,439,855,509,504
-30,545.54 -9,189.851 -355.7007 -24,769.72 -8,430.551 -73,291.37
-0.6338 -0.69147 -0.54455 -0.86282 -0.59012 -0.69714
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
366
M Santini, A Taramelli and A Sorichetta
considered PFTs. Though such differences are not higher than 1%, they are of the order of 104–106 hectares in absolute terms. Making a comparison between ASPHAA and the algorithm in LPJ at different resolutions and at global scale as in Section 3, their percent differences remain around 0.6%. An interesting conclusion is that the LPJ method appears to switch from over- to under-estimation when exceeding 15–20° latitude if compared with the SPHD and ASPH calculations. This reinforces the hypothesis about the usefulness of introducing latitude-based correction factors, especially to the most simplistic sphere-related approaches.
4.3.2 The CORINE land cover dataset The COordination foR Information oN the Environment (CORINE) project (http:// www.eea.europa.eu/publications/COR0-landcover) provided first in 1990 and then in 2000 the most detailed (at 1:100,000 scale) Pan-European land cover dataset, with accuracy higher than 85% (Feranec et al. 2010). Working on a smaller coverage with respect to global applications, the CORINE dataset for the year 2000 was here used to test the ASPHAA at finer resolutions. The dataset is available in a native vector format (http://www.eea.europa.eu/themes/landuse/clc-download) consisting of 100 km by 100 km tiles having a Lambert Azimuthal Equal Area (LAEA) PCS (see Appendix A3 for projection specification). All the tiles were first back-projected in GCS using the same Datum (ETRS1989) as the LAEA, and they were then mosaiced into four longitudinal zones as shown in Figure 9. Both projected and geographic versions of the vector maps for the four zones were rasterized, using a 500 m ¥ 500 m and a 0.005° ¥ 0.005° resolution for the former and the latter, respectively. While for the projected grids the total land use surfaces were computed considering a cell area equal to 250,000 m2, to the four lat/lon grids the ASPHAA was applied. Table 2 shows the results in terms of percent differences (diff %) and absolute percent differences (abs diff %) of aggregated land use areas calculated from PCS and GCS grids if compared with the reference surfaces calculated from source polygons. Either considering the whole European dataset or separately the four zones, the ASPHAA-based area has a better correspondence with the reference area, in terms of both relative and absolute differences. Results are influenced both by the accuracy of the source dataset and by the level of landscape fragmentation that causes biases during the vector to raster conversion, and for the ASPHAA-based areas also by the distance of the four zones from the origin of the projection of the source dataset. The mean centers (i.e. the geographic centers) for each land use polygon were first indentified, then their Euclidean distances from the center of LAEA projection (lon = 10°E and lat = 52°N) were computed and averaged when belonging to the same zone: as expected, in the case of absolute difference for ASPHAA, the biases increase when increasing the averaged mean center distance.
4.4 Regional Application The integration of grid-based biophysical parameters and high resolution demographic information at administrative level could be affected by the bias introduced through area © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
367
Figure 9 Map showing the four longitudinal zones in which the CORINE dataset was divided. The black point represents the center of the LAEA projection (lon = 10°E; lat = 52°N)
Table 2 Relative (diff %) and absolute (abs diff %) percent differences between total grid-based and vector-based (reference) surfaces Diff %
Abs diff %
Region
Projected
ASPHAA
Projected
ASPHAA
Averaged Mean Center distance from projection origin (m)
Whole dataset Zone 1 Zone 2 Zone 3 Zone 4
2.66 8.44 0.29 -1.63 -0.41
2.59 8.02 -0.03 -0.31 0.35
2.76 8.50 0.58 2.23 1.31
2.73 8.32 0.53 0.67 1.09
– 1,421,053 618,126 820,628 1,343,659
The averaged mean center distances represent the average distance of the mean centers for each land use polygon across the whole domain and within each zone from the origin of projection (lon = 10°E, lat = 52°N) © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
368
M Santini, A Taramelli and A Sorichetta
Table 3 Table shows the total area calculated with the two different coordinate systems (geographic and projected) within three classes of hazards related to Hurricane Stan Hazard class
Lat – Lon grid
UTM projection
1 2 3 Total affected area
24,053,670,000 pixels 8,026,730,000 pixels 5,687,519,000 pixels 37,767 km2
24,051,000,000 pixels 7,993,510,000 pixels 5,530,470,000 pixels 35,959 km2
calculation in a regular projected grid at the regional (sub-continental) scale. As an example of comparing and combining datasets from different sources and with different detail, a hurricane hazard/risk estimate in the Caribbean is here described for the Hurricane Stan event affecting these regions on September 2005 (Pasch and Roberts 2005). In this approach datasets about land use (http://bioval.jrc.ec.europa.eu/products/ glc2000/products.php), rain (http://www.cpc.noaa.gov/products/janowiak/cmorph_ description.html), wind speed (http://www.nhc.noaa.gov/pastall.shtml), topography (SRTM; http://seamless.usgs.gov/index.php) and bathymetry (GEBCO; http://www. gebco.net/) were processed and used for detecting areas vulnerable to storm surges, floods and high speed winds. Three multi-hazard severity classes (low, medium and high) were extracted by combining those representing the same severity classes referring to the single hazards, in order to assess the number of people affected by the considered event. Results are shown in Table 3, where is it possible to see the differences in number of pixels belonging to the three hazard severity classes, producing a total affected area (in km2) diverging when using geographic or projected (UTM, see Appendix A3 for projection parameters) coordinate systems, maintaining the source Datum (WGS 1984) and using 0.0008333333° ¥ 0.0008333333° and 90 m ¥ 90 m resolutions, respectively. The underestimation of the total affected area within the Caribbean by the projected coordinate system propagates the error to the final calculation of elements at risk in terms of affected population. Using population data at a 2.5′ ¥ 2.5′ angular resolution from the Gridded Population of the World dataset (GPWv3; http://sedac.ciesin.columbia.edu/ gpw/), the 1,418,634 affected people estimated using ASPHAA on the geographic grid is reduced to 1,010,210 affected people when projecting the grid using the UTM scheme for the zone 17, maintaining the source Datum (i.e. WGS 1984) and choosing a 5,000 ¥ 5,000 m linear mesh size spacing. This high underestimation by the projected grid could be related to the fact that UTM projection is a conformal and not an equal-area projection.
5 Conclusions With the advances of global position systems, remote sensing and GIS-computing technologies, there is an increasing need for practical, reliable and fast tools for spatial data analysis and processing. Either different spatial scales and resolutions, or biophysical and socio-economic data have to be considered in such evaluations. This allows a more in-depth analysis, through observation and modeling, of processes and phenomena occurring all over the globe, which are related to the dynamics of climate changes and economic development. © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
369
In this article, a well-developed algorithm, implemented into a GIS platform (ArcGIS9.x) and named ASPHAA, is proposed to calculate cell areas on gridded geodatasets in geographic coordinate systems and assuming a spherical shape for the Earth. Given several constraints due to byte representation for grid formats in ArcGIS and to the variety of resolutions that can be used under different research ambits, the analytical algorithm was first compared with a numeric solution of the cell area in the case of a more realistic spheroidal Earth’s model and considering seven different grid resolutions. Results confirm that: 1. the absolute percent difference does not exceed 0.01%; 2. the root mean square error is less than 30 ha in the worst case (coarser resolution, i.e. 10°); and 3. there is an increase of biases at higher latitudes. Then the algorithm was used at different spatial scales (from global to regional) and under different topics (land use assessment and scenarios, natural hazard risk evaluation, vegetation modeling) highlighting several points as listed below. Deviations from the true areas can be due first of all to vector-to-raster conversion performed when working at large scales to meet computational requirements. Moreover projected grids, both using equal-area (e.g. Goode Homolosine, Sinusoidal, LAEA) or other schemes (e.g. the conformal UTM) can create biases that, even if less than 2% and 4.5% at scales from global to regional, respectively, correspond to very large extents in absolute terms, ranging from millions of hectares to hundred thousands of hectares at scales from global to regional. The increase of percent deviations can be due to several reasons: aggregation from single to macro-categories in classifying data, the use of non-equal-area projections, the spread of the study area far from the projection point of the origin or the evaluations made by integrating data from different sources and with originally different spatial details. Also a recent study of Usery and Seong (2000) and Yildirim and Kaya (2008) showed for example how, using projected coordinate systems, total areas for several attributes is affected by projection type, latitude and with grid resolution; on the contrary very high performances of ASPHAA at any latitude and resolutions are here demonstrated. The comparison between ASPHAA and a simplified sphere-based algorithm without regard to authalic latitude conversion shows differences, whose sign is latitudedependent, which need to be considered for evaluating the uncertainty in applications where such algorithms are used. Although this study suggests that several improvements can be introduced in the presented algorithm, i.e. latitude- and resolution-based correction factors, others can be intuitively mentioned, as the accounting for a non-flat surface morphology of the cell, determining true areas certainly larger than from purely 2D estimations. At this stage, ASPHAA proved to be an interesting tool to improve multi-thematic and integrated research efforts through performing a more reliable area computation in multi-scale and multi-approach applications.
Appendix A1 – ASPHAA AML Code /*january 28th, 2009 /* command line: &r ASPHAA © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
370
M Santini, A Taramelli and A Sorichetta
&args grid &sv pi = 3.14159265358979 &describe %grid% &sv NROWS = %grd$nrows% &sv YLL = %grd$ymin%*%pi%/180 &sv CS = %grd$dx%*%pi%/180 &ty ‘Available spheroids:’ &ty ‘1 – AIRY 1830’ &ty ‘2 – AUSTRALIAN 1965’ &ty ‘3 – BESSEL 1841’ &ty ‘4 – BESSEL 1841 NAMIBIA’ &ty ‘5 – CLARKE 1866’ &ty ‘6 – CLARKE 1880’ &ty ‘7 – EVEREST 1830 INDIA’ &ty ‘8 – EVEREST 1830 MALAYSIA’ &ty ‘9 – EVEREST 1956 INDIA’ &ty ‘10 – EVEREST 1964 MALAYSIA & SINGAPORE’ &ty ‘11 – EVEREST 1969 MALAYSIA’ &ty ‘12 – EVEREST PAKISTAN’ &ty ‘13 – FISHER 1960’ &ty ‘14 – GRS 80’ &ty ‘15 – HELMERT 1906’ &ty ‘16 – HOUGH 1906’ &ty ‘17 – INDONESIAN 1974’ &ty ‘18 – INTERNATIONAL 1924’ &ty ‘19 – KRASSOVSKY 1940’ &ty ‘20 – MODIFIED AIRY’ &ty ‘21 – SOUTH AMERICAN 1969’ &ty ‘22 – WGS 72’ &ty ‘23 – WGS 84’ &ty ‘24 – ETRS 1989’ &sv spheroid:= [response ‘Choose the reference spheroid:’] &if %spheroid% = 1 &then &do &sv SMA = 6377563.396 &sv fl = 0.00334085067870327 &end &else &if %spheroid% = 2 &then &do &sv SMA = 6378160.0 &sv fl = 0.003352891899858333 &end &else &if %spheroid% = 3 &then &do &sv SMA = 6377397.155 &sv fl = 0.0033427731536659813 &end &else &if %spheroid% = 4 &then &do &sv SMA = 6377483.865 &sv fl = 0.003342773176894559 &end &else &if %spheroid% = 5 &then &do &sv SMA = 6378206.4 &sv fl = 0.0033900753039287908 &end
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm &else &if %spheroid% = 6 &then &do &sv SMA = 6378249.145 &sv fl = 0.003407561308111843 &end &else &if %spheroid% = 7 &then &do &sv SMA = 6377276.345 &sv fl = 0.0033244493186534527 &end &else &if %spheroid% = 8 &then &do &sv SMA = 6377298.556 &sv fl = 0.003324449343845343 &end &else &if %spheroid% = 9 &then &do &sv SMA = 6377301.243 &sv fl = 0.0033244493543833783 &end &else &if %spheroid% = 10 &then &do &sv SMA = 6377304.063 &sv fl = 0.003324449295589469 &end &else &if %spheroid% = 11 &then &do &sv SMA = 6377295.664 &sv fl = 0.0033244492833663726 &end &else &if %spheroid% = 12 &then &do &sv SMA = 6377309.613 &sv fl = 0.003324292418982392 &end &else &if %spheroid% = 13 &then &do &sv SMA = 6378155.0 &sv fl = 0.003352329944944847 &end &else &if %spheroid% = 14 &then &do &sv SMA = 6378137.0 &sv fl = 0.0033528106875095227 &end &else &if %spheroid% = 15 &then &do &sv SMA = 6378200.0 &sv fl = 0.0033523298109184524 &end &else &if %spheroid% = 16 &then &do &sv SMA = 6378270.0 &sv fl = 0.0033670034351006867 &end &else &if %spheroid% = 17 &then &do
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
371
372
M Santini, A Taramelli and A Sorichetta
&sv SMA = 6378160.0 &sv fl = 0.0033529256086395256 &end &else &if %spheroid% = 18 &then &do &sv SMA = 6378388.0 &sv fl = 0.003367003387062615 &end &else &if %spheroid% = 19 &then &do &sv SMA = 6378245.0 &sv fl = 0.0033523298336767685 &end &else &if %spheroid% = 20 &then &do &sv SMA = 6377340.189 &sv fl = 0.0033408506318589907 &end &else &if %spheroid% = 21 &then &do &sv SMA = 6378160.0 &sv fl = 0.003352891899858333 &end &else &if %spheroid% = 22 &then &do &sv SMA = 6378135.0 &sv fl = 0.00335277945669078 &end &else &if %spheroid% = 23 &then &do &sv SMA = 6378137.0 &sv fl = 0.0033528106718309896 &end &else &if %spheroid% = 24 &then &do &sv SMA = 6378137.0 &sv fl = 0.0033528106811823171 &end &sv e [sqrt [calc 2*%fl% - %fl%*%fl%]] &sv e1 [calc 1 - %e%*%e%] &sv e2 [calc 1/[calc 2*%e%]] &sv ee [calc %e%*%e%] &sv sin90 [sin [calc %pi%/2]] &sv m1 [calc 1 - %e%*%sin90%] &sv m2 [calc 1 + %e%*%sin90%] &sv mln [log [calc %m1%/%m2%]] &sv sin90_2 [calc %sin90%*%sin90%] &sv q [calc %e1%*[calc [calc %sin90%/[calc 1 - [calc %ee%*%sin90_2%]]] - [calc %e2%*%mln%]]] &sv R2 [calc %SMA%*%SMA%*[calc %q%/2]] &sv cost [calc %R2%*%CS%] grid setcell %grid% setwindow %grid% %grid% row = $$rowmap + 1 lgl = ((%NROWS% - row)*%CS%) + %YLL%/*lgl = lower geodetic latitude ugl = ((%NROWS% + 1 - row)*%CS%) + %YLL%/*ugl = upper geodetic latitude © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
373
q1 = %e1%*(sin(lgl)/(1 - %ee%*POW(sin(lgl),2)) - ((1/(2*%e%))*LN((1 %e%*sin(lgl))/(1 + %e%*sin(lgl))))) q2 = %e1%*(sin(ugl)/(1 - %ee%*POW(sin(ugl),2)) - ((1/(2*%e%))*LN((1 %e%*sin(ugl))/(1 + %e%*sin(ugl))))) cell_area = (%cost%*abs((q2/%q%) - (q1/%q%))) q &IF [exists row -grid] &THEN KILL row &IF [exists lgl -grid] &THEN KILL lgl &IF [exists ugl -grid] &THEN KILL ugl &IF [exists q1 -grid] &THEN KILL q1 &IF [exists q2 -grid] &THEN KILL q2
Appendix A2 – Numeric Solution for the Quadrangle Area on the Spheroid – Pseudo-code For each cell of the grid CA = 0 GLL = NROWS - NROW*CS GLLA = GLL GUL = (nrows-nrow+1)*CS GULA = GLL + 1/3,600 ACA = 0 For each arc-cell CLL = arctg (SMI2/SMA2)*tan(GLLA*pi/180.)) If it is not the highest latitude cell and arc-cell Then CUL = arctan (SMI2/SMA2)*tan(GULA* pi/180.)) Else CUL = pi/2 End If DCL = CUL - CLL RLL = sqrt((((SMA2*cos(GLLA*pi/180.))2) + ((SMI2*sin(GLLA*pi/180.))2))/ (((SMA*cos(GLLA*pi/180.))2) + ((SMI*dsin(GLLA*pi/180.))2))) RUL = sqrt((((SMA2*cos(GULA*pi/180.))2) + ((SMI2*sin(GULA*pi/180.))2))/ (((SMA*cos(GULA*pi/180.))2) + ((SMI*dsin(GULA*pi/180.))2))) !Using the Law of cosines (Carnot’s theorem) Dlat = sqrt(RLL2 + RUL2 - 2*RLL*RUL*cos(DCL) DLlow = RLL*cos(CLL)*((CS/nac)*pi/180.) DLupp = RUL*cos(CUL)*((CS/nac)*pi/180.) If it is not the highest latitude cell and arc-cell Then ACA= sqrt (Dlat 2-((DLupp - DLlow)/2)2)*(DLupp + DLlow)/2. Else ACA = sqrt (Dlat 2-((DLupp - DLlow)/2)2)*DLlow)/2. End If CA = CA+ACA GLLA = GLL + i* (1/3,600) If it is not the highest latitude cell and arc-cell Then GULA = 90. Else GULA = GLLA + (1/3,600) End If © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
374
M Santini, A Taramelli and A Sorichetta
End For CA = CA*NSL End For CA = Cell area GLL = Geodetic Lower Left cell corner latitude (in degrees) GLLA = Geodetic Lower Left corner latitude lower Arc-cell in the cell (in degrees) GUL = Geodetic Upper Left cell corner latitude (in degrees) GULA = Geodetic Upper Left corner latitude lower Arc-cell in the cell (in degrees) ACA = Arc-Cell Area (m2) CLL = geoCentric Lower Left corner latitude (in radians) CUL = geoCentric Upper Left corner latitude (in radians) DCL = difference between geocentric Upper and Lower Left corner latitude (in radians) RLL = Spheroid radius (m) at CLL RUL = Spheroid radius (m) at CUL Dlat = measure (in meters) of the meridian arc between upper and lower arc-cell corner DLlow = measure (in meters) of the parallel arc between left and right lower arc-cell corners DLupp = measure (in meters) of the parallel arc between left and right upper arc-cell corners Nac = number of arc-cells in a cell column i = order number of the arc-cells in the cell first column, counted from bottom to top NSL = number of columns with arc-cell horizontal resolution contained into a semilune with cell resolution.
Appendix A3 – Projection Parameters Projection: Goode_Homolosine False_Easting: 0.000000 False_Northing: 0.000000 Central_Meridian: 0.000000 Linear Unit: Meter (1.000000) Geographic Coordinate System: GCS_WGS_1984 Angular Unit: Degree (0.017453292519943299) Prime Meridian: Greenwich (0.000000000000000000) Datum: D_WGS_1984 Spheroid: WGS_1984 Semimajor Axis: 6378137.000000000000000000 Semiminor Axis: 6356752.314245179300000000 Inverse Flattening: 298.257223563000030000 Projection: Sinusoidal False_Easting: 0.000000 False_Northing: 0.000000 Central_Meridian: 0.000000 Linear Unit: Meter (1.000000) Geographic Coordinate System: GCS_WGS_1984 Angular Unit: Degree (0.017453292519943299) Prime Meridian: Greenwich (0.000000000000000000) Datum: D_WGS_1984 Spheroid: WGS_1984 Semimajor Axis: 6378137.000000000000000000 Semiminor Axis: 6356752.314245179300000000 © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
375
Inverse Flattening: 298.257223563000030000 Projection: Lambert_Azimuthal_Equal_Area False_Easting: 4321000.000000 False_Northing: 3210000.000000 Central_Meridian: 10.000000 Latitude_Of_Origin: 52.000000 Linear Unit: Meter (1.000000) Geographic Coordinate System: GCS_ETRS_1989 Angular Unit: Degree (0.017453292519943299) Prime Meridian: Greenwich (0.000000000000000000) Datum: D_ETRS_1989 Spheroid: GRS_1980 Semimajor Axis: 6378137.000000000000000000 Semiminor Axis: 6356752.314140356100000000 Inverse Flattening: 298.257222101000020000 Projection: Transverse_Mercator False_Easting: 500000.000000 False_Northing: 0.000000 Central_Meridian: -75.000000 Scale_Factor: 0.999600 Latitude_Of_Origin: 0.000000 Linear Unit: Meter (1.000000) Geographic Coordinate System: GCS_WGS_1984 Angular Unit: Degree (0.017453292519943295) Prime Meridian: Greenwich (0.000000000000000000) Datum: D_WGS_1984 Spheroid: WGS_1984 Semimajor Axis: 6378137.000000000000000000 Semiminor Axis: 6356752.314245179300000000 Inverse Flattening: 298.257223563000030000
Acknowledgements We thank Marco Traini for research assistance, and the referees for their critical reviews that helped to improve the manuscript. This work has been supported by the following grant: Fondazione Cassa di Risparmio di Foligno Award #1 Università di Perugia – Corso di Laurea in Protezione Civile. When this study began, Andrea Taramelli was affiliated with the LDEO of Columbia University and with the Department of Earth Sciences of the University of Perugia and Alessandro Sorichetta was affiliated with the Columbia University’s Center for International Earth Science Information Network (CIESIN).
References Albani M, Medvigy D, Hurtt G C, and Moorcroft P R 2006 The contributions of land-use change, CO2 fertilization, and climate variability to the Eastern US carbon sink. Global Change Biology 12: 2370–90 Arnold M, Chen R B, and Deichmann U 2006 Natural Disaster Hotspots Case Studies. Washington, DC, World Bank Publications Bach D, Barbour J, Macchiavello G, Martinelli M, Scalas P, Small C, Stark C, Taramelli A, Torriano L, and Weissel J 2007 Integration of the advanced remote sensing technologies to investigate the dust storm areas. In El-Beltagy A, Saxena M C, and Wang T (eds) Human and Nature: Working Tighter for Sustainable Development of Drylands. Aleppo, Syria, ICARDA: 387–97 © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
376
M Santini, A Taramelli and A Sorichetta
Bonham-Carter G F 1994 Geographic Information System for Geoscientist: Modelling with GIS. Oxford, Pergamon Burrough P A 1997 Environmental modeling with geographical information system. In Kemp Z (ed) Innovations in GIS. London, Taylor and Francis: 143–53 Burrough P A 2001 GIS and geostatistics: Essential partners for spatial analysis. Environmental and Ecological Statistics 8: 361–77 Crevoisier C, Shevliakova E, Gloor M, Wirth C, and Pacala S 2007 Drivers of fire in the boreal forests: Data constrained design of a prognostic model of burned area for use in dynamic global vegetation models. Journal of Geophysical Research 112: D24112, doi:10.1029/ 2006JD008372 Dilley M, Chen R B, Deichmann U, Lerner-Lam A L, and Arnold M 2005 Natural Disaster Hotspots: A Global Risk Analysis. Washington, DC, World Bank Publications Elvidge C D, Safran J, Tuttle B, Sutton P, Cinzano P, Pettit D, Arvesen J, and Small C 2007 Potential for global mapping of development via a nightsat mission. GeoJournal 69: 45–53 Engelstaedter S, Tegen I, and Washington R 2006 North African dust emissions and transport. Earth Science Reviews 79: 73–100 ESRI 1990 Understanding GIS: The ArcINFO Method. Redlands, CA, ESRI Press Feranec J, Jaffrain G, Soukup T, and Hazeu G 2010 Determining changes and flows in European landscapes 1990–2000 using CORINE land cover data. Applied Geography 30: 19–35 Food and Agriculture Organization (FAO) 2006a FAOSTAT 2005: FAO statistical databases [CD-ROM], Rome. WWW document, http://faostat.fao.org/default.aspx Food and Agriculture Organization (FAO) 2006b Agro-MAPS [CD-ROM], Rome. WWW document, http://www.fao.org/landandwater/agll/agromaps/interactive/page.jspx Fox-Rabinovitz M, Côté J, Dugas B, Déqué M, and McGregor J L 2006 Variable resolution general circulation models: Stretched-grid model intercomparison project (SGMIP). Journal of Geophysical Research 111: D16104, doi:10.1029/2005JD006520 Gringorten I I 1972 A square equal-area map of the world. Journal of Applied Meteorology 11: 763–67 Hay I, Graham A, and Rogers D (eds) 2006 Global Mapping of Infectious Diseases: Methods, Examples and Emerging Applications. London, Academic Press: 119–51 Henderson-Sellers A, Wilson M F, Thomas G, and Dickinson R E 1986 Current Global LandSurface Data Sets for Use in Climate Related Studies. Boulder, CO, National Center for Atmospheric Research Report No. NCAR/TN 272 Hertel T W 1997 Global Trade Analysis. Cambridge: Cambridge University Press Kimerling A J, Sahr K, White D, and Song L 1999 Comparing geometrical properties of global grid. Cartography and Geographic Information Science 26: 271–88 McCallum I, Obersteiner M, Nilsson S, and Shvidenko A 2006 A spatial comparison of four satellite derived 1 km global land cover datasets. International Journal of Applied Earth Observation and Geoinformation 8: 246–55 McGuire A, Sitch S, Clein J, Dargaville R, Esser G, Foley J, Heimann M, Joos F, Kaplan J, Kicklighter D, Meier R, Melillo J, Moore B, Prentice I, Ramankutty N, Reichenau T, Schloss A, Tian H, Williams L, and Wittenberg U 2001 Carbon balance of the terrestrial biosphere in the twentieth century: Analyses of CO2, climate and land use effects with four process-based ecosystem models. Global Biogeochemical Cycles 15: 183–206 Meshcheryakov G 1962 The best equal area projections. Geodesy and Aerophotogrametry 2: 132–7 Mitchell R B 2003 International environmental agreements: A survey of their features, formation and effects. Annual Review of Environment Resource 28: 429–61 Monfreda C, Ramankutty N, and Foley J A 2008 Farming the planet: 2, Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. Global Biogeochemical Cycles 22: GB1022, doi:10.1029/2007GB002947 Pasch R J and Roberts D P 2005 Tropical Cyclone Report – Hurricane Stan, 1–5 October 2005, National Hurricane Center. WWW document, http://www.nhc.noaa.gov/2005atlan.shtml Pasqui M, Tremback C J, Meneguzzo F, Giuliani G, and Gozzini B 2004 A soil moisture initialization method, based on antecedent precipitation approach, for regional atmospheric modeling system: A sensitivity study on precipitation and temperature. In Proceedings of the © 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)
ASPHAA: A GIS-Based Algorithm
377
Eighteenth Annual American Meteorological Society Conference on Hydrology, Seattle, Washington Pasqui M, Baldi M, Giuliani G, Gozzini B, Maracchi G, and Montagnani S 2005 Operational numerical weather prediction systems based on Linux cluster architectures. Nuovo Cimento Della Società Italiana Di Fisica C – Geophysics and Space Physics 28(2): 205–13 Perillo G M E, Piccolo M C, Mosquera J, and Aggio S 1999 Algorithm to calculate equal-area grid cells in irregular estuarine cross-sections. Computers and Geosciences 25: 277–82 Qihe Y, Snyder J P, and Tobler W 1999 Map Projection Transformation. New York, Taylor and Francis Ronneberger K, Tol R S J, and Schneider U A 2005 KLUM: A Simple Model of Global Agricultural Land Use as a Coupling Tool of Economy and Vegetation. Hamburg, Hamburg University and Centre for Marine and Atmospheric Science Ronneberger K, Berrittella M, Bosello F, and Tol R S J 2009 KLUM@GTAP: Introducing biophysical aspects of land-use decisions into a computable general equilibrium model a coupling experiment. Environmental Modeling and Assessment 14: 149–69 Sanyal J and Lu X X 2004 Application of remote sensing in flood management with special reference to monsoon Asia: A review. Natural Hazards 33: 283–301 Sitch S, Smith B, Prentice I C, Arneth A, Bondeau A, Cramer W, Kaplan J O, Levis S, Lucht W, Sykes M T, Thonicke K, and Venevsky S 2003 Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ Dynamic Global Vegetation Model. Global Change Biology 9: 161–85 Small C and Cohen J E 2004 Continental physiography, climate, and the global distribution of human population. Current Anthropology 45(2): 269–77 Snyder J P 1926 Map Projections: A Working Manual. Washington, DC, United States Geological Survey Professional Paper 1395 Thomas P 1952 Conformal Projections in Geodesy and Cartography. Washington, DC, United States Coast and Geodetic Survey Special Publication No 251 Tobler W R and Chen Z 1986 A quadtree for global information storage. Geographical Analysis 18: 4 Tobler W R, Deichmann R U, Gottsegen J, and Maloy K 1995 World population in a grid of spherical quadrilaterals. International Journal of Population Geography 3: 203–25 United Nations 1951 Annual Report on the International Map of the World of the Millionth Scale. New York, United Nations Social and Economic Affairs Division Usery E L and Seong J C 2000 A Comparison of Equal-Area Map Projections for Regional and Global Raster Data. In Proceedings of Twenty-ninth International Geographical Congress, Seoul, South Korea (available at http://carto-research.er.usgs.gov/projection/pdf/nmdrs.usery. prn.pdf) Washington R, Todd M C, Engelstaedter S, Mbainayel S, and Mitchell F 2006 Dust and the low-level circulation over the Bodélé Depression, Chad: Observations from BoDEx 2005. Journal of Geophysical Research 111: D03201, doi:10.1029/2005JD006502 Yetman G, Deichmann U, and Balk D 2000 Creating a Global Grid of Human Population. WWW document, http://sedac.ciesin.columbia.edu/gpw/ Yildirim F and Kaya A 2008 Selecting map projections in minimizing area distortions in GIS applications. Sensors 8: 7809–17 Yuan M 2005 Beyond mapping in GIS applications to environmental analysis. Bulletin American Meteorological Society 4: 169–70
© 2010 Blackwell Publishing Ltd Transactions in GIS, 2010, 14(3)