URBANIZATION AND ITS CARBON CONSEQUENCES IN THE YANGTZE RIVER DELTA, SOUTHEAST CHINA BY Ying Tang A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Geography 2011 ABSTRACT URBANIZATION AND ITS CARBON CONSEQUENCES IN THE YANGTZE RIVER DELTA, SOUTHEAST CHINA By Ying Tang China’s southeast coastal areas have undergone rapid expansion over the past few decades due to increasing population and economic growth. Although there are many studies focusing on the urban growth in this part of China, there is very little work that examines the carbon consequences of such growth. This research is concerned with the gross primary production change associated with urbanization process. In this study, MODIS Gross Primary Productivity (GPP) product, along with population, built-up area, and GDP data for 16 regions in the Yangtze River Delta were used to establish an empirical model which examined and quantified the relationship of urbanization and GPP change from 2003 to 2008. Also, Landsat remote sensing data calibrated with biophysical parameters were used to assess GPP using the light use efficiency (LUE) method for the Hangzhou city from 2001 to 2010. Additionally, the method of using vegetation phenology curves and climate indicators to forecast or backcast EVI was assessed. The empirical model results indicate that there is a significant positive association of a one-percent increase in built-up area ratio and a 15.4-percent decrease in annual total GPP for the Yangtze River Delta and capitals have a higher GPP given other things being equal. The primary cause of decreasing GPP is urban expansion, which encroached on large areas of surrounding cropland. Moreover, the estimation of GPP using phenology curves contains great uncertainty. The factors that strongly affect vegetation growth should be identified for a better estimation. To Mom and Dad iii Acknowledgements This work would not have been possible without support from many people. I owe my deepest gratitude to my advisor Dr. Jiaguo Qi, and my committee members, Dr. Ashton Shortridge and Dr. Nathan Moore, who provided thoughtful questions and suggestions. I am grateful to Dr. Lifeng Luo, Dr. Sue Grady, Dr. Kirk Goldsberry, Dr. David Lusch, Dr. David Campbell, Dr. Alan Arbogast, Dr. Bruce Pigozzi, Dr. Chen Guo, Dr. Shiyuan Zhong, Dr. Jessica Moy in MSU and Dr. Jiaping Wu in Zhejiang University. I’ve learned a lot from the courses they teach and benefited a lot from their encouragement and kind help. Thanks to the talk with Dr. Richard Groop and Dr. Jiaguo Qi during the “3+1” program interview in China. Otherwise, I would never choose geography as my major. Thanks to Shenpan Lin, Chuan Qin, Huiqing Huang, Miaoying Shi, Chenxiao Ling, Hui Xu, Xue Li, Siam Lawawirojwong, Michael Luehmann, and Paul McCord, who have been great colleagues and friends. Thanks to my boyfriend, my dad, my grandparents, and numerous friends, who offer love and support whenever I need them. Finally but not least, I want to thank my mom. My mom sacrificed everything for me. She worried about me all the time not because I don’t know how to take care of myself, but because no matter how old I grow, I’m always a child in her eyes. Without her support financially and emotionally, I would never be able to carry out this work. iv TABLE OF CONTENTS LIST OF TABLES ...................................................................................................................... vii LIST OF FIGURES ................................................................................................................... viii ABBREVIATIONS...................................................................................................................... xi 1 Introduction ............................................................................................................................ 1 1.1 Background ................................................................................................................. 2 1.2 Research Question ....................................................................................................... 5 1.3 Working Flow .............................................................................................................. 7 2 Literature Review .................................................................................................................. 8 2.1 Urbanization ................................................................................................................ 8 2.1.1 Conceptual Framework ..................................................................................... 8 2.1.2 Phases of Urbanization ..................................................................................... 8 2.1.3 Urban Remote Sensing ..................................................................................... 9 2.2 Gross Primary Production ......................................................................................... 12 2.2.1 Carbon Sequestration Measurement ............................................................... 12 2.2.2 Satellite-derived Measurement of Gross Primary Production ........................ 13 2.2.3 Impact of Urban Expansion on primary production ....................................... 14 2.3 Summary ................................................................................................................... 15 3 Study Area............................................................................................................................ 17 3.1 Yangtze River Delta .................................................................................................. 17 3.2 Hangzhou .................................................................................................................. 22 4 Data and Methods ................................................................................................................ 24 4.1 Changing GPP in Yangtze River Delta ...................................................................... 24 4.1.1 MODIS Land Cover Data ............................................................................... 24 4.1.2 MODIS GPP Product ...................................................................................... 26 4.1.3 Empirical Model ............................................................................................. 29 4.2 Changing GPP in the Hangzhou City ........................................................................ 32 4.2.1 Landsat-derived Land Cover Data .................................................................. 32 4.2.2 GPP derived using Light-use Efficiency Method ........................................... 35 4.3 Estimating GPP in a Data-sparse Environment using Phenology Curves ................. 40 4.3.1 Extracting Time-series Data and Seasonality ................................................. 40 4.3.2 Establish EVI-climate Function ......................................................................... 45 v 5 Results and Discussion ........................................................................................................ 48 5.1 Yangtze River Delta .................................................................................................. 48 5.1.1 Decreasing GPP in the Area ............................................................................ 48 5.2 Decreasing GPP in Hangzhou City ........................................................................... 58 5.2.1 Landsat-derived GPP ...................................................................................... 58 5.2.2 Comparison of Landsat GPP and MODIS GPP .............................................. 68 5.3 Method Derivation for Estimating EVI ..................................................................... 72 5.3.1 Vegetation Phenology in Hangzhou City ........................................................ 72 5.3.2 Function for Estimating EVI........................................................................... 76 6 Conclusion ........................................................................................................................... 84 Appendices.................................................................................................................................. 88 Bibliography ............................................................................................................................... 92 vi LIST OF TABLES Table 4.1 MOD 17 Product Output Layers ................................................................................ 27 Table 4.2 MODIS 17A2 Fill Values ........................................................................................... 27 Table 4.3 Summary Statistics of Variables for the Empirical Model ......................................... 32 Table 4.4 Landsat Image Used for Classification from 2001 to 2010 ....................................... 33 Table 4.5 BPLUT parameters for daily gross primary productivity .......................................... 39 Table 5.1 Summary Statistics of GPP in Yangtze River Delta ................................................... 55 Table 5.2 Summary of Regression Analysis for Equation (1) .................................................... 56 Table 5.3 Summary of Regression Analysis for Equation (2) .................................................... 57 Table 5.4 Summary Statistics of Landsat-derived GPP in Hangzhou City ................................ 61 Table 5.5 Summary of Temperature Change from 2001 to 2010 ............................................... 63 vii LIST OF FIGURES Figure 1.1 Working Flow of This Study ....................................................................................... 7 Figure 2.1 Cyclic model of the stages of urbanization based upon the population change in core and fringe zone of urban agglomerations: U, urbanization; S, suburbanization; D, disurbanization or counterurbanization; R, reurbanization phase (After Antrop, 2004) .............. 9 Figure 3.1 Map of Yangtze River Delta ...................................................................................... 18 Figure 3.2 Industry Percentage from 2003 to 2008 in the Yangtze River Delta ......................... 20 Figure 3.3 Population Change from 2003 to 2008 in the Yangtze River Delta .......................... 20 Figure 3.4 Built-up Area Change from 2003 to 2008 in the Yangtze River Delta ...................... 21 Figure 3.5 GDP Change from 2003 to 2008 in the Yangtze River Delta .................................... 21 Figure 3.6 The City of Hangzhou (Landsat Image Acquired in 2002) ....................................... 23 Figure 4.1 MODIS Sinusoidal Tiling System (after MODIS Land Team) ................................. 25 Figure 4.2 Landsat 7 SLC-off image (upper) and interpolated image (lower) ........................... 34 Figure 4.3 The TMIN and VPD attenuation scalars are simple linear ramp functions of daily TMIN and VPD........................................................................................................................... 39 Figure 4.4 Ten-year EVI Time-series Displayed in TIMESAT .................................................. 42 Figure 4.5 Ten-year EVI Time-series with Fitted Function ........................................................ 43 Figure 4.6 Processing Windows for Data-extracting .................................................................. 44 Figure 4.7 Seasonal Variation of EVI from 2001 to 2010 .......................................................... 46 Figure 5.1 GPP in Yangtze River Delta in 2003 ......................................................................... 49 Figure 5.2 GPP in Yangtze River Delta in 2005 ......................................................................... 50 Figure 5.3 GPP in Yangtze River Delta in 2007 ......................................................................... 51 Figure 5.4 GPP in Yangtze River Delta in 2009 ......................................................................... 52 viii Figure 5.5 Daily Average GPP in Yangtze River Delta from 2003 to 2010................................ 53 Figure 5.6 GPP Difference between 2003 and 2010 ................................................................... 54 Figure 5.7 Land Cover Type in Hangzhou in 2001 .................................................................... 59 Figure 5.8 Land Cover Type in Hangzhou in 2010 .................................................................... 60 Figure 5.9 Land Cover Change in Hangzhou from 2001 to 2010 .............................................. 60 Figure 5.10 GPP in Hangzhou City from 2001 to 2010 ............................................................. 62 Figure 5.11 Daily Mean Temperature from 2001 to 2010 .......................................................... 63 Figure 5.12 Daily T_scalar for Forest Land Cover from 2001 to 2010 ...................................... 64 Figure 5.13 Daily T_scalar for Crop Land Cover from 2001 to 2010 ........................................ 65 Figure 5.14 Daily VPD_scalar for Forest Land Cover from 2001 to 2010 ................................ 66 Figure 5.15 Daily VPD_scalar for Crop Land Cover from 2001 to 2010 .................................. 67 Figure 5.16 Mean Daily GPP Derived from Two Algorithms in Hangzhou City ....................... 71 Figure 5.17 The relationship of Landsat GPP and MOD 17 GPP .............................................. 72 Figure 5.18 Time-series EVI for Built-up Area .......................................................................... 73 Figure 5.19 Time-series EVI for Forest Land Cover Type ......................................................... 74 Figure 5.20 Time-series for Crop Land Cover Type ................................................................... 75 Figure 5.21 The Relationship of EVI and Temperature (T0 to T3) ............................................ 77 Figure 5.22 The Relationship of EVI and Temperature (T4 to T7) ............................................ 78 Figure 5.23 The Relationship of EVI and Precipitation (P1 to P4) ............................................ 79 Figure 5.24 The Relationship of MODIS EVI and Precipitation (P5 to P8) .............................. 80 Figure 5.25 The Relationship of Function-fitted EVI and Temperature (T1 to T2) ................... 81 ix Figure 5.26 Estimated EVI from 2001 to 2010 .......................................................................... 81 x ABBREVIATIONS APAR Absorbed Photosynthetically Active Radiation AVHRR Advanced Very High Resolution Radiometer BPLUT Biome Parameter Look-Up Table DMSA Detroit-Ann Arbor-Flint Metropolitan statistical Area ETM+ Enhanced Thematic Mapper Plus EVI Enhanced Vegetation Index FPAR Fraction of Photosynthetically Active Radiation GDP Gross Domestic Product GLDAS Global Land Data Assimilation System GPP Gross Primary Production IGBP International Geosphere Biosphere Programme LAI Leaf Area Index LUE Ligh Use Efficiency LULC Land Use Land Cover LUT Look Up Table MODIS MODerate resolution Imaging Spectroradiometer NASA National Aeronautic and Space Administration NDVI Normalized Difference Vegetation Index NIR Near-Infrared xi NOAA National Oceanic and Atmospheric Administration NPP Net Primary Production OLS Ordinary Least-Squares SLC Scan Line Corrector TM Thematic Mapper VI Vegetation Index VPD Vapor Pressure Deficit xii 1 Introduction During the past sixty years, humankind has experienced dramatic population growth from 2.5 billion to 6.7 billion. About sixty-three percent of this gain took place in urban areas, especially in developing countries. Over the past three decades, China has experienced rapid development from an economic standpoint. With more than 19 percent of the world’s population and the world’s second largest Gross Domestic Product (GDP), China’s land use and land cover (LULC) change will certainly have effects across national borders and deserves global attention. Human activities change the landscape to a great extent which causes a significant influence on carbon cycling. This research aims to quantitatively examine the relationship between urbanization and Gross Primary Production (GPP) change in one of China’s most economically vigorous regions – the Yangtze River Delta. Moreover, different approaches of estimating GPP are compared and a new method of using phenology curves to estimate GPP in a data-sparse environment is established. For this part of study, the city of Hangzhou was chosen as an example study site. The objectives of this study are threefold: 1) establish an empirical model to assess the relationship between GPP change and urbanization in the Yangtze River Delta; 2) derive higher-resolution GPP in the Hangzhou city using light-use efficiency model, and compare derived GPP to MODIS GPP; 3) establish the relationship between Vegetation Index (VI) and climate indicators to facilitate VI estimation in a data-sparse environment. The remainder of Chapter 1 discusses the background of this study and its main purpose. Chapter 2 concentrates on the review of the development of studies on urbanization and GPP, as 1 well as their relationship. Chapter 3 provides a brief introduction of the study region: the Yangtze River Delta and Hangzhou City. Chapter 4 has three subsections. The first section focuses on urbanization and its carbon consequences for the entire Yangtze River Delta. An empirical model is constructed with socioeconomic data and GPP value derived from MODerate resolution Imaging Spectroradiometer (MODIS) GPP product. The second section aims to look at GPP and LULC change in the Hangzhou city. GPP is derived using Light-Use Efficiency (LUE) method with Landsat images and MODIS Enhanced Vegetation Index (EVI) products. Two algorithms for GPP calculation are compared and analyzed. For the third section, the main purpose is to develop and implement an algorithm to forecast or backcast EVI in the Yangtze River Delta based on vegetation phenology and climate indicators from 2001 to 2010. The Hangzhou City again is used as an example. The research results and discussion are presented in Chapter 5. Chapter 6 summarizes this research and discusses directions extending from the presented study. 1.1 Background When the People’s Republic of China was formally established in 1949, there were only 132 1 cities in the country. This number increased to 655 by the end of 2008, indicating China’s rapid urbanization rate. The urban population reached 0.62 billion with a national urbanization rate of 46.6 percent by 2009. From 2000 to 2009, urban population increased by 0.163 billion. A report released by the Chinese Academy of Social Sciences in 2010 predicted that the urban population will reach 52% in 2015 and swell to 65% by 2030. In 2010, China’s GDP reached 5,878.6 billion dollars and replaced Japan as the world’s second-largest economy. Along with such growth, 1 Based on reports from National Bureau of Statistics of China. 2 China’s land use and land cover have also undergone dramatic change, southeastern China in particular. Urbanization has great impacts on carbon cycling. In developing countries, economic growth attracts young people who live in rural areas to quit agricultural activities and move to cities. Rapid population growth requires new residential areas. With the development of industry, more commercial space is needed. As a result of these processes, cities expand their boundaries by encroaching upon surrounding lands occupied by agricultural activities and natural vegetation. On the other hand, economic growth may generate greater demand in luxury goods such as meat (Sicular, 1985), which requires more land to engage in agricultural activities or more intense use of existing agricultural land. Urbanization also has great influence on regional land-atmosphere interaction, energy balance, and regional climate. Urban heat island effect changes the surrounding area’s temperature and atmosphere (Bounoua et al., 2009). Humans alter atmospheric composition by burning fossil fuel. Urban impervious materials cannot hold water like soil; they alter the water cycle greatly. All those processes have direct or indirect influences on the carbon cycle. The most basic needs of food, fiber, and fuel which humankind survive upon are all supplied by biological productivity. Consequently, the dynamics of carbon cycling of terrestrial ecosystem call for great attention globally. Urbanization process can affect carbon cycling positively or negatively. By seeking the linkage between urbanization and biological productivity change, insights will be provided for planning and natural resource management for both short-term and long-term. 3 The rate at which plants convert light energy during photosynthesis process is termed primary production. There are two types of primary production used to describe the amount of the fixed carbon by vegetation. Gross primary production is the sum of all the energy converted, which indicates the uptake of carbon from the atmosphere; net primary production (NPP) equals GPP less the energy loss during plant respiration, which is associated with the net flux of carbon between vegetation and atmosphere (Zhao, 2007). Remote sensing has been widely used to estimate GPP and NPP due to its ability to detect and measure photosynthetic energy, vegetation volume as represented by vegetation indices, and light conversion efficiency factor of different types of vegetation (Running, 2004). GPP’s measurement is more directly related to remote sensing measurements and thus is chosen in this study to quantify primary production change. Spatial and temporal resolution and extents are important components of this work. For a regional study of GPP and LULC change, MODIS products can be adopted. For studies focusing on individual cities, the coarse resolution of MODIS products may introduce substantial errors and uncertainty. In these cases, Landsat images are more suitable to be used to derive LULC maps. For calculation of GPP after 2000, the MODIS VI products enable data inputs with a high temporal resolution. However, before 2000, the MODIS VI product is not available and only one or two cloud-free Landsat images are usable per year for classification purposes in this region of China. To more accurately assess GPP change, more extensive temporal coverage is needed. To address this issue, I derive phenology curves for different vegetation types and establish the relationship between VI and precipitation, temperature with the facilitation of phenology 4 parameters in order to model GPP based on seasonal biological characteristics of plants. The possibility of adopting this approach is examined. 1.2 Research Question This study considers three research questions. First, this study aims to establish the relationship of urbanization and GPP change in the Yangtze River Delta from 2003 to 2008. To achieve that, urbanization process needs to be quantified using several indicators (i.e. regional GDP, buit-up area, population, etc.); summation of GPP will be calculated for each region every year; and finally, the sum of GPP will be regressed against the urbanization indicators. Second, as MODIS GPP product is designed to accommodate large geographic area analysis, studies for individual cities should adopt different approaches. In this part, a Landsat-based light-use efficiency model is empolyed to derive GPP for the Hangzhou city. GPP for Hangzhou is derived with MODIS EVI, Landsat images, ground-station climate data and surface shortwave incident radiation data. This modeled GPP will be referred to as Landsat GPP in the following chapters. Its algorithm has MODIS EVI as input, and the term Landsat GPP is chosen only to distinguish it from the MODIS GPP product. Landsat GPP outputs will be compared with the MODIS 17A2 GPP product. Third, the potential to use EVI phenology curve to model GPP in a data-sparse environment (i.e. one usable Landsat image per year) is assessed. Phenology curves are a useful tool to assess vegetation growth. By associating this curve to gross primary production, it is possible to assess GPP in China from 1970s to 2000, prior to the launch of the MODIS sensor and using only one or two cloud-free Landsat images per year. Great insights may be provided toward carbon 5 accounting and long-term policy making. For this part, the analysis will also be performed using the Hangzhou city as an example. 6 Socioeconomic Data 1.3 Working Flow Built-up Area Assess Urbanization’s Carbon Consequences in the Yangtze River Delta Empirical Model Population GDP Enable Small Scale GPP Analysis with Finer Resolution Landsat Image MODIS GPP Atmospheric Correction Reprojection Stripe Interpolation for SLC-off Image Unsupervised Classification MODIS EVI Landsat Image BPLUT GPP in Hangzhou City LULC FPAR GLDA Data Climate Data LUE Method Shortwave Radiation Tscalar Min T Max T Dew Point T Temperature Epsilon es VPD ea Regression Analysis EVI Estimation F(t) Facilitate GPP Estimation in a Data-sparse Environment Figure 1.1 Working Flow of This Study 7 2 Literature Review 2.1 Urbanization 2.1.1 Conceptual Framework Although there is currently no international definition for urbanization, the most widely accepted definition express urbanization as the process of the concentration of population to cities or suburbs. From the prospective of economy, cities are the centers where people engage in non-agricultural activities. Thus urbanization is viewed as the transformation process from agricultural activity to non-agricultural activity. From the angle of anthropology, urbanization is a process which changes human life style. For geographers, besides the change of population and economy, urbanization is also concerned as a spatial process. Generally, urbanization can be reflected in the following aspects: (1) concentration of population, including both expansion of scale and increase of concentration centers; (2) increase of proportion of urban population in total population; (3) prevalence of urban life style and (4) change from agricultural activity to non-agricultural activity. 2.1.2 Phases of Urbanization Being a long temporal evolution process, urbanization has several phases: urbanization, suburbanization, counterurbanization and reurbanization (Champion, 2001). Those four stages correspond to stages of population migration from fringe to center, faster suburb growth than inner city growth, decline of population in both center and fringe and recovery of population in city center and then in fringe, accordingly (Antrop, 2004) (Fig. 1). 8 Figure 2.1 Cyclic model of the stages of urbanization based upon the population change in core and fringe zone of urban agglomerations: U, urbanization; S, suburbanization; D, disurbanization or counterurbanization; R, reurbanization phase (After Antrop, 2004) Over the past decade, China has experienced rapid urbanization while its largest cities have rapidly suburbanized (Cervero and Day, 2008). Ever since 1950s, many developed countries started to experience suburbanization. In China, suburbanization is relatively new. Nanjing, Shanghai, and Hangzhou are among the first batch of cities that experienced suburbanization in the Yangtze River Delta. The suburbanization process has placed many residents in locations far away from the city centers (Ma, 2004). Economic growth, transportation facility development, urban planning and policy making are shown to be the key drivers of suburbanization (Jing, 2007). 2.1.3 Urban Remote Sensing During the past three decades, remotely sensed data have been applied to numerous urban studies (Welch, 1980; Haack et al., 1986; Gallo et al., 1993; Jensen et al., 1999; Weng, 2002; Pan et al., 2007). Although maps remain a valuable source of information, remote sensing products have advantages in providing multi-spatial resolution and multi-temporal resolution information 9 to support different objectives of various researches. Moreover, for less developed areas where project funds are not abundant, publicly accessible remote sensing products (Landsat, MODIS, etc) become major data sources. Applying remote sensing to urban areas is relatively new compared to work focused on the natural environment (Weng et al., 2006). With the development of technology and social needs, urban remote sensing has received more and more interest. To meet the special characteristics of urban landscapes, new methods, algorithms, technologies, and products were developed. To extract urban features, LIDAR and SAR could be utilized to generate large-scale urban orthoimage (Rottensteiner et al, 2002; Dell et al., 2001). Urban landscapes are typically composed of features smaller than a pixel from medium resolution satellite imagery, such as Landsat Thematic Mapper (TM) or Enhanced Thematic Mapper Plus (ETM+). Subpixel analysis algorithms and bayesian spectral mixture analysis thus could be adopted under this situation to produce fraction image with a more realistic representation of the nature of a surface (Ji et al, 1999; Song, 2005). Besides urban land use and land cover monitoring, detection and prediction, remote sensing is also used in applications of planning and socioeconomic studies including urban heat island identification, environmental assessment, and population estimation (Streutker, 2003; Güneralp et al., 2008; ). In China, the first city to conduct an aerial remote sensing experiment was Tianjing (Chen et al., 2000). National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) and Landsat images were used to produce ecological environment atlas of Beijing-Tianjin-Tangshan region (Fu et al., 1989). During 1985 to 1995, a series of studies were carried out with remote sensing images on environmental issues over more 10 than 90 cities in China (Chen et al., 2000). With thirty years’ efforts, remote sensing has contributed greatly to land resource planning, decision-making and management of various projects of the State Land Management Bureau, such as nation-wide second land use investigation, the state key mineral resources investigation, and ecological environment investigation in the Qinghai-Tibet Plateau, the Yangtze River Valley, and the Yellow River valley. Coastal cities that are medium to large in size in southeastern China and large cities in the inner land draw most attention in urban expansion researches (Seto et al., 2005; Xiao et al., 2006; Weng, 2002). Weng combined remote sensing, GIS and stochastic modeling in assessing land use and land cover changes in the Zhujiang Delta, China. During the period 1989 – 1997, urban area and horticulture farms increased in area, with a decrease in area of cropland (Weng, 2002). Seto et al. examined urban growth for Guangzhou, Zhoushan, Donguan and Shenzhen using time series landscape metrics derived from remote sensing images. Results show that the average annual rate of urban land-use change was 17% between 1988 and 1999. Total urban land for the four cities almost quadrupled during the study period (Seto et al., 2005). Xiao et al. investigated urban expansion and land use change in Shijiazhuang with remote sensing and GIS. The urban 2 2 area of Shijiazhuang city expanded from 6.31 km in 1934 to 165.5 km in 2001 at an average 2 rate of 2.4 km per year (Xiao et al., 2006). The Yangtze River Delta has also experienced great land use change during the past three decades: main expansion areas are widely distributed along and to the south of Yangtze River, presenting the trend of urban agglomerations; Shanghai, Nanjing, Suzhou-Wuxi-Changzhou, Hangzhou-Shaoxing-Ningbo, Zhenjiang-Yangzhou-Taizhou formed as five centers of land use 11 dynamic areas (He et al., 2006). He also concluded that large cities basically expanded from the developed center; medium cities developed along main transportation routes; cities that are relatively small developed along the main transportation lines toward the regional developing center. Urbanized areas were mostly converted from farmland, followed by conversion of a portion of forestland and grassland (Seto et al., 2005; He et al, 2006; Weng, 2002; Seto et al., 2003; Xu et al., 2000). Population growth, GDP growth, and the State's macro-control and economic policy play significant roles in driving urbanization. 2.2 Gross Primary Production 2.2.1 Carbon Sequestration Measurement Various methods have been established to estimate carbon sequestration. Biomass field measurement, eddy-covariance flux tower method, model estimation, and remote sensing method are most common ones found in literature (Niu et al., 2008). Field measurement is the most direct method. It can be conducted by tree diameter measurement, growth rate analysis, defoliation collection and so forth. The objective of field measurement is to calculate above-ground net primary productivity and underground net primary productivity (Ni, 2004). However, its high cost in time and labor limits the applications towards large study areas. For eddy-covariance method, CO2, H2O, wind and temperature are measured to quantify atmospheric flux and to determine ecosystem gross primary production, autotrophic respiration, heterotrophic respiration and net ecosystem exchange (Curtis et al., 2002). Ecosystem models simulate the carbon cycle of energy, nutrients and micronutrients in order to provide estimations of carbon 12 budget. Even though the modeling process compensate for the discreteness of measurement, accuracy of modeling is still bounded by ecological mechanisms at multi-scales. Date collected from remote sensing satellites contribute to the possibility of estimating carbon balance over large areas with intensive temporal repetition. By transforming from multi-points data to continuous surface data, spatial distribution pattern over the study area could also be examined. 2.2.2 Satellite-derived Measurement of Gross Primary Production Before the 1980s, biologists focused on the organism level and ecological studies, and no attention was paid to global scale study (Running et al., 2004). The earliest attempt was made by geographers to investigate ecological process at global scale. In 1975, Lieth and Whittaker published first global NPP estimates. These estimates were based on regressions of actual evapotranspiration calculated from temperature against a few NPP field plots (Lieth and Whittaker, 1975). Running et al. (2004) concluded that three activities built the foundation of global scale terrestrial ecology: the expansion of atmospheric flask sampling usage enabled gaining information on a biospheric scale about the photosynthetic uptake and evolution of CO2 (Tans et al., 1990); more submodels of land surface processes were incorporated into global climate models; efforts were made by a few ecologists to apply meteorological satellite data on vegetation analysis. The first image of Normalized Difference Vegetation Index (NDVI) became the cover of Science in August 1985 which also gave impetus to global vegetation research (Tucker et al., 1985). 13 Monteith was the first to propose the relationship between Absorbed Photosynthetically Active Radiation (APAR) and NPP with the logic that NPP of well-watered and fertilized annual crop plants was linearly related to APAR (Monteith, 1972; Monteith, 1977). In the early 1980s, National Aeronautic and Space Administration (NASA) was encouraged by the Reagan administration to pursue ecological science projects with a broader geographic scope. With such impetus, new exploration trying to relate NDVI to GPP, NPP began to emerge. Goward et al. (1985) examined seasonal variation of NDVI and proved its agreement with vegetation phenology. Running et al. analyzed growing-season dynamics and related NDVI to simulation of photosynthesis and transpiration (Running et al., 1988). These studies presented the great potential in applying remote sensing to ecological process. The light use efficiency method was established based on previous theories and experiments and widely accepted and adopted in numerous studies (Turner et al., 2003; Xiao et al., 2005; Stavros et al., 2007; Wang et al., 2010). MODIS Terra was launched in 1999 and the MODIS Aqua satellite was launched in 2002. Both satellites are operating across the world every day and distributing GPP products every 8 days globally. The theoretical basis of light use efficiency method and algorithm for global GPP will be discussed in the Chapter 4. 2.2.3 Impact of Urban Expansion on primary production Vegetation productivity is the source of all food, fiber and fuel for humankind consumption. Rapid urban sprawl and associated dense population and economic development have exerted high pressure on ecosystem. 14 In the United States, urban land transformation in the mid-1990s has reduced the amount of carbon fixed through photosynthesis by 0.04 petagram per year (Imhoff et al., 2004). Their results also show that even though urbanization is taking place on the most fertile land, the reduction in NPP due to urbanization of agricultural lands is equivalent to the caloric requirement of 16.5 million people. The change of GPP associated with urbanization has notable spatial difference. The Detroit-Ann Arbor-Flint Metropolitan statistical Area (DMSA) of southeastern Michigan has an increased GPP during 1991 and 1999. It is mainly due to the increase in the fraction of tree cover throughout the entire region (Zhao et al., 2007). Lu et al. examined the effects of population, GDP, and settlement on GPP in southeastern China with Pearson’s correction analysis based on remote sensing data in 2000. The analysis shows that NPP is negatively correlated with settlement, population and GDP with settlement generally linearly related to NPP and population and GDP nonlinear related to NPP (Lu et al., 2010). In Shenzhen, 80.1% change of land use is due to the urban sprawl and cropland reduction during 1999 – 2005 with a decrease of total NPP from 1811.0 Gg C to 1489.49 Gg C (Yu et al., 2009). Other studies in China also show similar results that urbanization leads to reduction in GPP or NPP (Xu et al., 2007; Hu et al., 2009). 2.3 Summary In consideration of large scale analysis on the effects of urbanization on GPP, adopting remote sensing data is necessary and advantageous. Current research on subjects like urbanization or NPP, GPP change in China usually provides single snapshots based on one image or examine two to three time points; regression analysis are also always based on these one to 15 three time points; most research calculated NPP or GPP based on NDVI or Leaf Area Index (LAI) (Imhoff et al., 2004; Yu et al., 2009; Xu et al., 2007). With the MODIS GPP product, seasonal variability of GPP could be examined. Moreover, MODIS GPP 8-day composite output can be used to show seasonal variability in analysis, which may provide more insights than simply use one summed value for each year, providing more samples in the dataset as well. As MODIS GPP product is not designed for small scale analysis, GPP calculation based on Landsat image may provide important information on GPP analysis for individual cities. A method combining a recent algorithm which adopts Enhanced Vegetation Index as model input will be compared with original MODIS algorithm. More details about this method will be discussed in the method chapter. 16 3 Study Area 3.1 Yangtze River Delta The Yangtze River Delta, when referred to geographically, is the area that was formed from gathered mud and sand brought down by the river. It is generally thought to be the territory east of Zhenjiang city, south of Tongyang Canal, and north of Hangzhou Bay. It is a part of Middle-lower Yangtze Plain. As now people are more concerned about the socioeconomical status of this area, the geographical boundary of the Yangtze River Delta has three different definitions. The first definition is referred to as Small Yangtze River Delta which includes 16 regions including Shanghai, southern Jiangsu province and northern Zhejiang province. These 16 regions are Shanghai, Wuxi, Ningbo, Zhoushan, Suzhou, Yangzhou, Hangzhou, Shaoxing, Nanjing, Nantong, Changzhou, Huzhou, Jiaxing, Zhenjiang, Taizhou (pingyin, Tàizhou, referred to as Taizhou4 in later descriptions to distinguish with the other city) and Taizhou (pingyin, Tāizhou, referred to as Taizhou1 in the following chapters). In 1992, 14 regions that are close in geographic locations organized the meeting in which the 14-region cooperation committee was founded. In 1997, with the joining of Taizhou4 into this association, the name of the committee was changed to Yangtze River Delta Urban Economic Coordination Committee. The event that Taizhou1 joined the committee in 2003 is a milestone that the “Yangtze River Delta” represents an economic zone more than only a geographic area. With 6 other regions applied to join the committee, the committee has 22 region members by 2010. The State Council issued “On The Yangtze River Delta Region To Further Promote The 17 Reform And Opening Up And Economic And Social Development Of Guidance” in 2008 which brought in the second definition of Yangtze River Delta – all the area of Shanghai, Jiangsu Province and Zhejiang Province. The third definition introduces the idea that Anhui Province and Jiangxi Province should also be included in the Yangtze Economic Zone. Figure 3.1 Map of Yangtze River Delta (For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis.) In this study, the Yangtze River Delta only refers to the 16 regions mentioned in the first definition (Figure 3.1). It is mainly because of their close geographic locations and similar 18 economic status. Also, the study period spans mainly before 2010. So regions joined the committee in 2010 are not considered in this study. Data assemblage issues in different versions of statistical yearbooks is another reason that only the 16 regions are chosen as study region. Another important point to get here is that the 16 regions that was referred to earlier include 16 cities and their subordinate cities and counties, while further analysis of Hangzhou in the following chapters only refer to districts of Hangzhou City, not including its subordinate cities and counties. In other words, whenever Hangzhou is mentioned together with other 15 regions, it is referring to Hangzhou City and its subordinate cities and counties, as shown in Figure 3.1; otherwise, it only means the city itself (Figure 3.6). The Yangtze River Delta Region exhibits a “W” shape with Shanghai at the joint of two “V”s, cities from the Jiangsu Province at the north wing and cities from the Zhejiang Province at the south wing. It covers a large area of 109,890 km2. In 2009, this region contributed to 20% of national GDP with approximately 1% of total land area and 7% of total population of China. The GDP contribution is much greater than Pearl River Delta and Jing-Jin-Tang area, showing the leading position of Yangtze River Delta in economic growth. During the recent years, the Yangtze River Delta is experiencing a decline in the primary industry, an increase in tertiary industry, and the secondary industry is relatively stable (Figure 3.2). The State government encourages greater development of tertiary industry and further adjustment of industry structure. 19 0.6 Industry Percentage 0.5 0.4 0.3 Primary Industry 0.2 Secondary Industry 0.1 Tertiary Industry 0 2003 2004 2005 2006 2007 2008 Year Figure 3.2 Industry Percentage from 2003 to 2008 in the Yangtze River Delta Population (10,000 people) 8450 8400 8350 8300 8250 8200 8150 8100 2003 2004 2005 2006 2007 2008 Year Figure 3.3 Population Change from 2003 to 2008 in the Yangtze River Delta 20 4000 Built-up Area (km2) 3500 3000 2500 2000 1500 1000 500 0 2003 2004 2005 2006 2007 2008 Year Figure 3.4 Built-up Area Change from 2003 to 2008 in the Yangtze River Delta 6000 GDP (billion Yuan) 5000 4000 3000 2000 1000 0 2003 2004 2005 2006 2007 2008 Year Figure 3.5 GDP Change from 2003 to 2008 in the Yangtze River Delta The Yangtze River Delta has a marine monsoon subtropical climate. The typical weather is hot and humid in summer, cool and dry in winter, and warm in spring and fall. This region contains the most fertile soil in China. While paddy rice is the dominant crop in this region, there are also large areas of cropland growing cotton, wheat, peanut and rapeseed. 21 The Yangtze River Delta is one of the regions have best natural resources in China. However, the rapid economic growth has high pressure on the environment and natural resources, which threats sustainable development of this area. Among all the issues, land use is one primary issue which is often neglected during development. From 1980 to 1995, the Yangtze River Delta lost 2 247,000 hm agricultural lands with a losing rate 6.7 times higher than the State average. With the loss of agricultural lands, the remained agricultural lands may suffer more intense use with stronger use of fertilizer, which may introduce degradation issues in the future. The growing population, fast economic development and industry structure adjustment will all have great influences on LULC change and ecosystem (Figure 3.3-3.5). This makes Yangtze River Delta a great region to study urbanization and its carbon consequences. Furthermore, to study the method of derivation of GPP for a smaller geographic area, the city of Hangzhou was chosen as an example. 3.2 Hangzhou Hangzhou is the capital city of Zhejiang Province. It has a long history of more than 2200 years since the Qin Dynasty governed it as a county. It’s located in the south wing of Yangtze River Delta, west of Hangzhou Bay. The center geographic coordinates of the city is 30°16´N, 2 120°12´ E. There are 8 districts in the city, with a total area of 3068 km (Figure 3.6). According to the census in 2008, there are 7.97 million people living in the city, of which 69.3% is urban population. The population density is 480 persons per square kilometer. Comparing to 2007, the ratio of Hangzhou’s three industries in 2008 changed from 4.1:50.2:45.7 to 3.7:50.0:46.3. The center of the city is built-up area, surrounded by large area of cropland. Forests are distributed 22 near the western boundary. The Hangzhou City is chosen due to its rapid development and urbanization, various land cover types (crop and forest), and more available data. Figure 3.6 The City of Hangzhou (Landsat Image Acquired in 2002) 23 4 Data and Methods 4.1 Changing GPP in Yangtze River Delta This section will discuss the preparation process of building an empirical model to address the relationship of urbanization and GPP: calculating GPP data for each one of 16 region in the Yangtze River Delta, gathering and choosing socioeconomic data factors that have influence on GPP change, and construct the equation to describe the relationship of two subjects. 4.1.1 MODIS Land Cover Data MODIS is a key instrument aboard the Terra and Aqua satellites. Terra and Aqua are two satellites viewing the earth every one to two days with Terra passing the equator from north to south in the morning and Aqua passing the equator from south to north in the afternoon. They acquire data in 36 spectral bands at three spatial resolutions: 250m, 500m and 1000m. The MCD12Q1 yearly product was used in this study. It is derived from observations of both Terra and Aqua with 500m grid. The MODIS Land Cover Type product contains five classification schemes. The primary classification scheme used in this study identifies 17 land cover types which are defined by the International Geosphere Biosphere Programme (IGBP). 24 Figure 4.1 MODIS Sinusoidal Tiling System (after MODIS Land Team) 25 The MODIS Land Cover Type product uses the Sinusoidal grid tiling system with tiles that are 10 degrees by 10 degrees in size at the equator. The tile coordinate system starts at (0, 0) in the upper left and ends at (35, 17) in the lower right corner (Figure 4.1). To encompass the whole area of Yangtze River Delta, four tiles of images are needed. They are h27v5, h27v6, h28v5 and h28v6. To enable processing the files in ERDAS, ENVI and ArcGIS, MODIS Reprojection Tool was used to reproject and mosaic these images. Data from 2001 to 2009 were downloaded, reprojected to geographic coordinate system WGS1984 and reclassified into five types: built-up, crop, tree, water and others. 4.1.2 MODIS GPP Product The MODIS GPP product is a cumulative estimation of GPP based on radiation use efficiency. Each output is an 8-day composite at 1 kilometer spatial resolution. There are 46 th output images per year. The 46 output is the summation for 5 days in a non-leap year or 6 days in a leap year. The accuracy of GPP products is assessed and the uncertainties are analyzed taking global conditions into consideration. The Psn_QC_1km layer contained in the MODIS GPP product is produced to assure data quality (Table 4.1). However, in this study, no GPP data was ruled out, due to the consistency in the control of regression analysis. One assumption is that all cities have somewhat similar chance of getting good quality GPP estimates due to their close geographic locations. 26 Table 4.1 MOD 17 Product Output Layers Science Data Sets (HDF Layers) Gpp_1km: Gross Primary Production Units PsnNet_1km: Net Photosynthesis Kg PSN_QC_1km: QC for GPP/PSN Kg Bit Type Fill Valid Range 2 16-bit signed 32761–32767 integer 2 16-bit signed 32761–32767 -30000–30000 integer C/m C/m Bit 8-bit unsigned integer 255 0–30000 Multiply By Scale Factor 0.0001 0.0001 0–254 NA This study used MODIS 17A2 to calculate GPP for the Yangtze River Delta area. The MOD 17 products are also in the integerized sinusoidal projection. Data for four tiles (h27v5, h27v6, h28v5 and h28v6) were downloaded, reprojected, and mosaiced at 1 km resolution. As most of standard MODIS products have a valid range, the valid data range for MODIS 17 GPP layer is 0 – 30000. Any value beyond the range is fill value representing different types of land cover (Table 4.2). Batch command was executed in ERDAS for all files to change values that fall out of this range to zero. All the images were kept in integer format with values 10000 larger than actual index value. Table 4.2 MODIS 17A2 Fill Values Value 32767 32766 32765 32764 32763 32762 32761 Description Fill value land cover assigned as perennial salt or inland fresh water land cover assigned as barren, sparse vegetation (rock, tundra, desert.) land cover assigned as perennial snow, ice land cover assigned as "permanent" wetlands/inundated marshlands land cover assigned as urban/built-up land cover assigned as "unclassified" or not able to determine 27 The processed GPP product was then subset to the extent of 16 regions, i.e. Shanghai, Wuxi, Ningbo, Zhoushan, Suzhou, Yangzhou, Hangzhou, Shaoxing, Nanjing, Nantong, Changzhou, Huzhou, Jiaxing, Zhenjiang, Taizhou4 and Taizhou1. The area with fill values were replaced with data derived using MODIS algorithm with NDVI-FPAR function2. FPAR stands for the Fraction of Photosynthetically Active Radiation (FPAR). Two values were summarized for each 8-day composite image for each region. The first one is the total sum of GPP calculated by adding raster values for all cells together, which will be referred to as GPP sum in the following discussion. GPP sum quantifies the amount of GPP produced in certain cities each year. The second value is the mean GPP value of all the cells. Equation for converting GPP sum into GPP mean can be written as: ∑ġ GPPsum(i) † GPP sum = GPP mean * Cell(number) Annual GPP sum = Annual GPP density = Annual GPP Sum / Cell (number) 2 GPP mean: 0.0001 kg C/(m * year) Cell (number): the total number of raster cells that fall in the city boundary The cell number is not the product of column number and row number because of the irregular shapes of cities. If the product of column number and row number was used, GPP mean would be underestimated. In order to count the cell numbers for each city, raster layers were converted to ASCII format 2 Method can be found on Algorithm Theoretical Basis Document , http://modis.gsfc.nasa.gov/data/atbd/atbd_mod15.pdf 28 and summarized. It is also important to note that the GPP sum value is not in the standard unit. In order to get the total GPP in unit of kg C for each study area, the size that each cell represents 2 (i.e., 1 km ) should be multiplied to GPP sum. 4.1.3 Empirical Model The empirical work considers the impact of urban expansion on GPP through a balanced panel of 16×6=96 observations (one for each city-year from 2003 to 2008 within the 16 cities located in Yangtze River Delta). The two subsections give the analysis of effect on GPP from two perspectives: first build an annual GPP sum equation, and then construct an annual GPP mean equation. These two functions can estimate the average effect of urban expansion on the quantity of GPP in this area. The built-up area and built-up ratio are used as two urbanization indicators in the equations mentioned above. They measure the expansion process directly. Population and GDP are also considered in the empirical model: a larger population often means higher pressure on the ecosystem (i.e. disturbance of natural environment and vegetation growth); regions with higher GDP may be more willing to sacrifice environment protection for economic growth or more focused on “smart growth”, which preserves natural environment and make planning decision more wisely. 4.1.3.1 Annual GPP Sum Equation This subsection investigates the effect of urban expansion on annual GPP sum through equation (1): (1) log GPPIT = α0 +a1builtup IT +hM' IT + lT' + α2 state2 + α3 state3+ ε1 29 where GPP IT represents the annual GPP sum for city I and year T; ε1is the random error items for equation (1). The dependent variable, logarithm of annual GPP sum in the function, depends on: 1) builtup IT is the urbanization indicator for city I and T; 2) M' IT vector other characteristics including logarithm of population, logarithm of land area, logarithm of GDP, and the dummy variable of whether it is a capital city; 3) T' vector 5 year dummy variables for each year from 2004 to 2008; 4) state 2 and state 3 are the province dummy variable for Jiangsu and Shanghai, respectively. 4.1.3.2 Annual GPP Density Equation This subsection considers the effect of urbanization on annual GPP density through equation (2). (2) log GPPDENSITYIT = b 0 + b1builtup IT + uN' IT + dT' + b 2 state2 + b 3 state3 + ε 2 where GPPDENSITY IT represents the annual GPP density for city I and year T; ε2 is random error items for equation (2). The dependent variables, logarithm of annual GPP density in the function, depends on: 1) builtup IT is the urbanization indicator for city I and T; 2) N' IT vector other characteristics including logarithm of population density, logarithm of GDP, and the dummy variable of whether it is a capital city; 3) T' vector 5 year dummy variables for each year from 2004 to 2008; 4) state 2 and state 3 are the province dummy variable for Jiangsu and Shanghai, respectively. 4.1.3.3 Socioeconomic Data and Urbanization Indicators This subsection describes the data that is used to perform the analysis in equation (1) and (2). As discussed above, the key independent variable is the urban expansion indicator. The 30 econometric model including equation (1) and (2) presented above is examined two times, once for each urban expansion indicator. The first eligible indicator, builtup - 1 IT , measures the logarithm of built-up area for city I and year T. The second indicator, builtup - 2 IT , measures the percentage of built-up area in total land area. As a result, the coefficient of the first indicator can predict the effect when built-up 2 2 area increased by one percent (i.e. the built-up area increased from 100 km to 101 km ) and the coefficient of the second indicator can predicts the effect when built-up area ratio increased by 2 one percent (i.e. the land area for a city is 1000 km , within which the built-up area ratio 2 2 increased from 10 percent (100 km ) to 11 percent (110 km ). MODIS LULC product has a fixed urban/built-up boundary and is thus not suitable for urbanization analysis. So efforts were made to collect indicator parameters from statistical reports. Both the data of built-up area (BUILTUPAREA) and the land area (LANDAREA) for each city-year are available in THE YEARBOOK OF CHINA’S CITIES. Therefore, the mathematical formula of builtup - 1 IT and builtup - 2 IT are: builtup - 1 IT = log( BUILTUPARE AIT ) builtup - 2 IT = BUILTUPARE AIT ´ 100 LANDAREA IT THE YEARBOOK OF CHINA’S CITIES also provides the data the following variables including city population, and GDP that play as independent variables in equation (1) and (2). The population density can be derived through dividing population by land area in a city. Summary statistics of the data used in the regression analysis are provided in Table 4.3. 31 Table 4.3 Summary Statistics of Variables for the Empirical Model Variable GPP Mean 4.30×10 Std. Dev. 7 7 3.64×10 Min Max 6 7.70×10 1.79×10 8 Variable Explanation and Unit Cell number * 0.0001 Kg 2 C/(m * year) GPPDENSITY 6552.2 2120.0 3570.1 11768.9 BUILTUPAREA 187.229 204.162 34 886 LANDAREA 6864.95 3319.46 1440 16596 Builtup-1 Builtup-2 POPULATION POPDENSITY GDP capital 4.83638 0.822889 3.52636 6.78672 2.93652 3.087871 0.712411 13.97476 518.1245 275.9652 96.58 1391.04 0.0801702 0.038677 0.038731 0.2194069 2.37×107 2.40×107 1.7×106 1.40×108 0.1875 0.3923613 0 1 state2 0.4375 0.4986825 0 1 state3 0.0625 0.2433321 0 1 2 0.0001 Kg C/(m * year) 2 Km 2 Km percentage 10,000 people 10,000 people / Km2 10,000 RMB =1 if it is a capital city; otherwise=0 =1 if it is in Jiangsu Province; otherwise=0 =1 if it is Shanghai; otherwise=0 4.2 Changing GPP in the Hangzhou City 4.2.1 Landsat-derived Land Cover Data Landsat MSS/TM/ETM+ data collected spanning 1970s to 2010 were downloaded from NASA website. The entire Hangzhou city is located within path 119 and row 39 of Landsat swath. The satellite images of this region often have large area of cloud cover. Due to this reason, it is usually the case that only one scene is qualified for use of classification each year. Even though it is favorable to have images acquired at approximately the same time of every year, such requirement is often impossible to achieve. From 2001 to 2010, one image per year was 32 classified to get the land cover information (Table 4.4). During this study period, most of the images were acquired by Landsat ETM+ with one exception – classification of 2005 used a mosaic image from Landsat TM and ETM+ due to cloud coverage. Table 4.4 Landsat Image Used for Classification from 2001 to 2010 Acquired Time 2001/02/16 2002/07/13 2003/03/26 2004/02/09 2005/10/17, 2005/11/26 2006/04/19 2007/05/08 2008/01/03 2009/03/10 2010/12/10 Sensor Landsat ETM+ Landsat ETM+ Landsat ETM+ Landsat ETM+ Landsat TM and ETM+ Landsat ETM+ Landsat ETM+ Landsat ETM+ Landsat ETM+ Landsat ETM+ As Landsat level 1T data, each image has been geographically registered in WGS 1984 UTM projection UTM Zone 50 with a cell size of 30 meter. Each band was then atmospherically corrected using the method in the paper written by Chander et al., which summarized radiometric calibration coefficients for Landsat MSS, TM and ETM+. On May 31, 2003, an instrument malfunction occurred onboard Landsat 7. The problem was caused by failure of the Scan Line Corrector (SLC), which results in the stripes in SLC-off Landsat data. Many efforts and progresses have been made to assess the usability and develop 3 interpolating method for Landsat 7 product ever since . In this study, Landsat images with a 3 USGS SLC-off usability assessments, 2003 33 SLC-off issue were interpolated with an add-on function of ENVI software. The Landsat SLC-off image and the interpolated result are shown in Figure 4.2. As suggested by previous research, even though interpolated images may not be as accurate as the SLC-on image, they remain valuable in providing LULC information and calculating statistics, which provides more useful information than masked images. After that, all the bands except the thermal bands were stacked together and clipped to the extent of Hangzhou city for subsequent analysis. Figure 4.2 Landsat 7 SLC-off image (upper) and interpolated image (lower) 34 Unsupervised classification was performed to classify images into five categories: built-up, crop, tree, water, and others. The raster images were first grouped into 80 classes and classified based on false color composite images. The first round of classification usually leaves many mixed features. Then, those pixels could not be classified in the first round will be exported for unsupervised classification again. It usually took six to ten cycles to classify a single image. During the classification process, Google Earth was used to verify as the ground truth as it has a finer resolution and also has historical images at multiple time points. The atmospheric correction process and unsupervised classification process were both conducted using ERDAS IMAGINE software. Following classification, the urban boundary was delineated combining classification results and visual interpolation with false color image on the background. 4.2.2 GPP derived using Light-use Efficiency Method 4.2.2.1 Light-use Efficiency Method Satellite-based studies have been using light-use efficiency approach to estimate GPP (Running et al., 2000; Imhoff et al., 2004; Stagakis et al., 2007). Significant efforts and progress have been made during the past decade. The method used in this study is a combined approach based on the MODIS GPP algorithm with some adjustments from a recent algorithm published by Xiao et al.. Running et al. proposed the GPP algorithm based on Monteith’s logic that the primary production of well-watered and fertilized annual crop plants is related to the amount of absorbed solar energy linearly. So the daily GPP could be calculated by: GPP = ε * APAR 35 APAR = PAR * FPAR ε=εmax * TMIN_scalar * VPD_scalar where APAR is absorbed photosynthetically solar energy and ε is the energy-to-carbon conversion efficiency (i.e., LUE, g C / MJ), PAR is the daily incident radiation, and FPAR is the fraction of PAR absorbed by vegetation. TMIN and VPD scalars are temperature and Vapor Pressure Deficit (VPD) linear attenuation scalars that quantify the effect of severe climatic constraints. Subfreezing temperature and high vapor pressure deficits may force leaf stomata to close. εmax can be found from Biome Parameter Look-Up Table (BPLUT). In the original algorithm developed by Running et al., FPAR is equal to or approximately equal to normalized difference vegetation index (NDVI) which uses red band and near-infrared (NIR) band reflectance (Running et al., 2000). NDVI = (ρnir - ρred)/ (ρnir + ρred) EVI is calculated as: EVI = G * (ρnir - ρred) / (ρnir + C1* ρred – C2 * ρblue + L) where G = 2.5, C1=6, C2=7.5, and L=1. The ρblue, ρred, and ρnir represent reflectance at the blue (0.45-0.52µm), red (0.6-0.7µm), and Near-Infrared (NIR) wavelengths (0.7-1.1µm), respectively (Huete et al., 1997). MODIS GPP product adopts MOD 15 algorithm for FPAR, which is also based on NDVI-FPAR relationship. Grigera and Oesterheld carried out a study in wheat and pasture sites to compare parameterization of FPAR and LUE. Their study suggested that even though that NDVI-FPAR relationship was adopted in many satellite products’ algorithms and widely used, 36 enhanced vegetation index (EVI) may be more suitable for deriving FPAR as it presented a more linear relationship with FPAR at study sites (Grigera and Oesterheld, 2006). Xiao et al. used EVI in their studies on forests and grasslands and suggested a better correlation of EVI to GPP than NDVI (xiao et al., 2005). In their study, FPAR is calculated as: FPAR = a * EVI a=1 Comparing to NDVI, EVI algorithm corrects some distortions in the reflected light caused by the particles in the air and the ground cover below the vegetation. It also does not become saturated as easily as the NDVI in densely vegetated areas like rain forests. This study combines light-use efficiency method from previous studies and calculates GPP by: GPP =εmax * TMIN_scalar * VPD_scalar * EVI * PAR PAR = SWRad * 0.45 TMIN_scalar = (TMIN – TMIN_min) / (TMIN_max – TMIN_min) VPD_scalar = (VPD_max – VPD) / (VPD_max – VPD_min) where SWRad is incident shortwave radiation, TMIN is daily minimum temperature, VPD can be calculated using daily average temperature and relative humidity. Parameters for TMIN_scalar and VPD_scalar are illustrated in Table 4.5 and Figure 4.3. VPD is the difference of mean saturation vapour pressure (es) and actual vapour pressure (ea), es and ea can be calculated with daily maximum temperature, daily minimum temperature, and dew point temperature (Td): 37 VPD = es – ea ÚúǴƼ es(T) = 0.6108 * es = Ǭ exp a ∗Ú⁄(Ú Ǭxp. ) .x (Úú† ) ea = es(Td) = 0.6108 * exp a ∗Ú ⁄(Ú Ǭxp. ) .x In this study, TMIN_scalar and VPD_scalar are calculated with equations stated above based on ground-station data of daily maximum temperature, daily minimum temperature, daily average temperature, and daily dew-point temperature from the NOAA National Climatic Data Center. The ground station is located at N 30°13′58″E 120°10′1″, in the Shangcheng District of Hangzhou City. Generally, spatial interpolation is favorable in estimating temperature surface. However, stations in surrounding cities are geographically too far away and some stations only have data before 2000 (i.e. Changzhou, etc). So data from station No. 584570 was used in calculation. This approach may introduce uncertainty in ε estimation. The results of daily TMIN_scalar and VPD_scalar were aggregated into 23 outputs for each one-year period. The ε value for built-up area was assessed based on averaged built-up green percentage (36%) assuming an even distribution of trees and grasses in the built-up land and estimated percentage of 18% and 18%. The ε value of tree, crop/grass and built-up types for each one of 23 periods are stored in integer format multiplied by scale factor of 1000,000. 38 Table 4.5 BPLUT parameters for daily gross primary productivity Parameter εmax TMINmax Units (kgC/MJ) (°C) TMINmin (°C) VPDmax (Pa) VPDmin (Pa) Description The maximum radiation conversion efficiency The daily minimum temperature at which ε = εmax (for optimal VPD) The daily minimum temperature at which ε = 0.0 (at any VPD) The daylight average vapor pressure deficit at which ε = 0.0 (at any TMIN) The daylight average vapor pressure deficit at which ε = εmax (for optimal TMIN) Figure 4.3 The TMIN and VPD attenuation scalars are simple linear ramp functions of daily TMIN and VPD In this study, EVI is from MODIS 13 vegetation indices product and has a 250m spatial resolution and a 16-day temporal resolution. For MODIS Terra, EVI product has 23 outputs per year starting from Day 001 (i.e., 001, 017, 033, …, 353). The εmax for forest and crop are from 4 the Biome Properties Look-Up Table (BPLUT) . SWRad data are from Global Land Data Assimilation System (GLDAS) Noah Land Surface 4 MODIS 17 User Guide, http://www.ntsg.umt.edu/modis/MOD17UsersGuide.pdf 39 Model L4 monthly 0.25 x 0.25 degree product in NetCDF format. They are downloaded, converted to raster, and aggregated into 16-day composite by assigning certain number of days from a certain month to the intended 16-day period output (i.e., 001, 017, … , 353). 4.3 Estimating GPP in a Data-sparse Environment using Phenology Curves With the facilitation of various remote sensing products, estimating GPP remotely for a large area has never been easier. However, most of these products are enabled after 2000. The physical world has changed a lot since the 20 th century. Past always provides significant meaning to the future. It would be of great value to be able to assess GPP in the 1970s. This section will introduce the data and method used to derive vegetation phenology curve and the possibility of adopting this approach in carbon storage assessment. 4.3.1 Extracting Time-series Data and Seasonality MODIS 16-day 250m global VI product (i.e. MOD 13Q1) is used in this study to extract time-series and seasonality parameters. There are ten-year (2001 to 2010) EVI available with 23 images per year. TIMESAT software developed by Lars Eklundh and Per Jönsson was used in data processing. Raster images of EVI were converted to flat binary file with IDL programming. Converted binary files were then viewed, processed in TIMESAT. Outputs containing time-series were written to ASCII files for further processing (Figure 4.4). To reduce noise (cloud contamination effect, etc.) in the data, original time-series data were fitted with double logistic functions. The fitted data were processed and written to ASCII files afterwards (Figure 4.5). The 40 seasonal parameters were extracted to better understand the vegetation phenology in Hangzhou City. The seasonality parameters include the start of the season, the end of the season, and the length of the season. The EVI image for Hangzhou city has 333 rows and 393 columns within which only selected windows (each window contains a few rows and columns) with forest land cover were examined to develop the algorithm. The processing windows are chosen by editing column and row numbers in TIMESAT (Figure 4.6). The seasonality parameters were extracted within the whole image extent. 41 EVI Number of 16-day Periods Figure 4.4 Ten-year EVI Time-series Displayed in TIMESAT 42 EVI Number of 16-day Periods Figure 4.5 Ten-year EVI Time-series with Fitted Function 43 Figure 4.6 Processing Windows for Data-extracting 44 4.3.2 Establish EVI-climate Function Iwasaki tried to predict NDVI distribution over Mongolia with precipitation and temperature data a few months before prediction (Iwasaki, 2009). The grids were aggregated into spatial resolution of 0.25°×0.25°and analyzed with stepwise linear regression. The predicted results agree well with the actual NDVI. This method may well suit the need for short-term prediction. However, long-term forecast or backcast may bring in more issues such as human disturbance, ecosystem evolution, etc. So this study will only perform the analysis on the relative stable land cover type, i.e. the forest land cover. To estimate EVI for a relatively long time period, the de-trended value for each cell’s EVI is first presented with a function (Figure 4.7): F(t) = C1cos(wt) + C2sin(wt) + C3 where w = 20π/t, and t = 1, 2, 3, …, 230. EVI_diff = EVI – F(t) EVI_diff2 = EVI_fit – F(t) where EVI_diff is the difference between MODIS EVI and F(t) function. 45 0.8 0.6 0.4 0.2 0 1 23 45 67 89 111 133 155 177 199 221 Figure 4.7 Seasonal Variation of EVI from 2001 to 2010 The next step is to establish the relationship of the difference between actual EVI and F(t) and climate conditions. As F(t) function represents the seasonal cycle of vegetation growth, the main purpose is to test whether it is possible to explain interannual variation through temperature or precipitation difference. There are some assumptions of made in adopting this method: the land cover type did not change significantly during the study period; the vegetation species remains stable during the study period; seasonal differences are mainly attributing to climate conditions. As vegetation growth is affected by temperature and precipitation a few months before, Iwasaki established NDVI function for one grid as: NDVI_Aug = 0.04 ×Precpitation_Jul + 0.002 × Temperature_Mar + 1.01 Efforts are made in this study to build the function in the format of: EVI_diff = B1×Precpitation(tn) + B2 × Temperature(tm) + B3 EVI_diff2 = B4×Precpitation(tn) + B5 × Temperature(tm) + B6 46 where EVI_diff and EVI_diff2 are differences between actual EVI or fitted EVI with F(t) function, B1, B2, B3, B4, B5, and B6 are coefficients, tn and tm are time periods before EVI acquisition. Each cell would have a set of coefficients across the ten-year period. 47 5 Results and Discussion 5.1 Yangtze River Delta 5.1.1 Decreasing GPP in the Area GPP for each one-year period in Yangtze River Delta were added together and summarized (Table 5.1). GPP sum for several years are shown in Figure 5.1 - 5.4. Figure 5.5 presents the 10-year 8-day daily average GPP in Yangtze River Delta Region. Each 8-day period GPP was divided by 8 with 5-day or 6-day periods divided by 5 or 6. The X axis shows time increment 2 and the Y axis is GPP mean in units of 0.0001 Kg C/(m *day). The higher values correspond to the time of growing season of vegetation while the lowest GPP values are often found during the winter. From 2003 to 2010, the decreasing GPP trend is obvious. Figure 5.6 shows the difference of annual GPP of 2010 and 2003. The northern area of Yangtze River Delta presents a higher GPP of 2010 than 2003 while most of the southern part has a lower GPP. For the whole delta, the mean GPP in the nonurban area decreased by 12.5 percent from 2003 to 2010. The average annual decrease is approximately 1.8 percent. The empirical model in the next section presents the quantitative relationship and GPP change and urbanization. 48 0 Figure 5.1 GPP in Yangtze River Delta in 2003 49 50 100 Km 0 Figure 5.2 GPP in Yangtze River Delta in 2005 50 50 100 Km 0 Figure 5.3 GPP in Yangtze River Delta in 2007 51 50 100 Km 0 Figure 5.4 GPP in Yangtze River Delta in 2009 52 50 100 Km 50 45 Daily mean GPP 40 35 30 25 20 15 10 5 0 1 177 353 526 702 875 1051 1225 1401 1574 1750 1923 2099 2273 Day Figure 5.5 Daily Average GPP in Yangtze River Delta from 2003 to 2010 53 2449 2622 2798 High: 7840 0 Low: -12147 0 50 Figure 5.6 GPP Difference between 2003 and 2010 54 100 Km Table 5.1 Summary Statistics of GPP in Yangtze River Delta Year 2003 2004 2005 2006 2007 2008 2009 2010 Mean 6709.464 6832.934 6100.812 6164.190 6365.581 5987.022 6193.568 5867.574 Min 0 0 0 0 0 0 0 0 Max 27641 26643 27926 24955 25339 25096 25401 22244 Std. Dev. 4820.187 4466.003 3913.752 4053.594 4378.622 4074.264 4410.703 4107.373 2 Mean, Min, Max: 0.0001 kg C/(m *year) 5.1.2 Empirical Model Results Table 5.2 and Table 5.3 report Ordinary Least-Squares (OLS) regressions of equation (1) and (2) using the balanced panel of 16×6=96 observation, respectively. 5.1.1.1 Results from Equation (1) Column 1 in Table 5.2 shows that 1) a one-percent increase in built-up area (i.e. the built-up area increased from 100 km2 to 101 km2) is significantly correlated to 1.066 percent decrease in annual GPP sum; 2) a one-percent increase in population is significantly correlated to a 1.176-percent decrease in annual GPP sum; 3) a one-percent increase in land area is strongly correlated to a 1.368-percent increase in annual GPP sum; 4) a one-percent increase in GDP is strongly correlated to a 0.994-percent increase in annual GPP sum; and 5) capital city has a 1.468 times higher annual GPP sum on average compared with other cities when other things being equal. Column 2 in Table 5.2 shows that 1) a one-percentage increase in built-up area ratio (i.e. the land area for a city is 1000 km2, within which the built-up area ratio increased from 10 percent (100 km2) to 11 percent(110 km2) ) is significantly correlated to a 15.4-percent decrease in annual GPP sum; 2) a one-percent increase in land area is strongly correlated to a 0.865-percent increase in annual GPP sum; and 3) capital city has a 0.902 times higher annual GPP sum on average compared with other cities when other things being equal. 55 Table 5.2 Summary of Regression Analysis for Equation (1) Dependent Variable log (GPP) (1) Urbanization Indicator Builtup-1 -1.066*** (0.198) Builtup-2 Other Independent Variable log (POPULATION) log (LANDAREA) log (GDP) capital Year fixed effects year2004 year2005 year2006 year2007 year2008 (2) -0.154** (0.054) -1.176*** (0.308) 1.368*** (0.246) 0.994*** (0.213) 1.468*** (0.263) -0.544 (0.330) 0.865* (0.343) 0.273 (0.154) 0.902** (0.268) -0.0846 (0.154) -0.216 (0.158) -0.390* (0.172) -0.489* (0.186) -0.612** (0.200) 0.0504 (0.168) -0.0459 (0.170) -0.0923 (0.178) -0.114 (0.187) -0.171 (0.197) Province fixed effects state 2 -0.162 (0.132) state 3 0.0475 (0.365) _cons 1.297 (1.936) R-square 0.6959 No. of Obs. 96 Standard Error in parentheses, * p < 0.05, ** p < 0.01, *** p < 0.001 -0.116 (0.145) 0.672 (0.425) 8.916*** (1.601) 0.6262 96 5.1.1.2 Results from Equation (2) Column 1 in Table 5.3 shows that 1) a one-percent increase in built-up area (i.e. the built-up area increased from 100 km2 to 101 km2) is significantly correlated to a 56 0.155-percent decrease in annual GPP density; 2) a one-percent increase in population density is strongly correlated to a 0.248-percent decrease in annual GPP density; and 3) capital city has a 20.0 percent higher GPP density on average compared with other cities when other things being equal. Column 2 in Table 5.3 shows that 1) a one-percentage increase in built-up area ratio (i.e. the land area for a city is 1000 km2, within which the built-up area ratio increased from 10 percent (100 km2) to 11 percent (110 km2)) is significantly correlated to a 4.92-percent decrease in annual GPP density; and 2) capital city has a 22.2 percent higher GPP density on average compared with other cities when other things being equal. Table 5.3 Summary of Regression Analysis for Equation (2) Dependent Variable log (GPPDENSITY) (1) Urbanization Indicator Builtup-1 -0.155* (0.063) Builtup-2 Other Independent Variable log (POPDENSITY) log (GDP) capital Year fixed effects year2004 year2005 year2006 year2007 year2008 Province fixed effects state 2 (2) -0.0492*** (0.013) -0.248* (0.094) 0.0755 (0.046) 0.200* (0.088) -0.0682 (0.105) -0.0221 (0.028) 0.222** (0.074) 0.0446 (0.057) -0.0365 (0.058) -0.0649 (0.060) -0.0656 (0.062) -0.109 (0.064) 0.0713 (0.055) -0.00126 (0.063) -0.0112 (0.056) 0.0022 (0.055) -0.0289 (0.060) 0.409*** (0.049) 0.423*** (0.048) 57 Table 5.3 (cont’d) state 3 0.122 (0.137) _cons 7.400*** (0.587) R-square 0.7622 No. of Obs. 96 Standard Error in parentheses, * p < 0.05, ** p < 0.01, *** p < 0.001 0.288* (0.135) 8.823*** (0.524) 0.7806 96 In summary, the increase of build-up area and population is strongly related to GPP decrease: the increasing built-up area may take the space of cropland and forest, causing overall less vegetation in the whole area; on the other hand, the increasing needs for food along with the population growth is not satisfied or supplied by the Yangtze River Delta, food are importing from elsewhere, so cropland area didn’t have large increase; higher population density may lead to the degradation of ecosystem in this area. Capitals like Nanjing, Shanghai, Hangzhou have higher GPP mean than other cities given equal population and GDP. It is probably because the capital cities focus more on the protection of environment and turn away high potential pollution industries. The capitals are usually key tourism cities: they are usually more concerned about increasing green space and attracting more people for eco-tourism. Given other things being equal, the positive correlation of GDP and GPP may because that regions with a better economic status are more concerned about environment protection, or simply have more funds for planting trees and protecting forests. Two equations actually established the relationship of urbanization and GPP change from four models. Even though the regression results do not have identical parameters, the general information that they carry are similar. 5.2 Decreasing GPP in Hangzhou City 5.2.1 Landsat-derived GPP Hangzhou experienced rapid urban growth during 2001 to 2010. While forest land cover 58 remained stable, large area of crop were converted to built-up land (Figure 5.7-5.9). GPP for Hangzhou city was calculated based on light-use efficiency method for 2001 – 2010 period with MODIS 16-day EVI product, Landsat image classification results, and surface shortwave incident radiation data from GLDA. Each 16-day GPP (13-day for the last cycle of non-leap year, 14-day for leap year) was then accumulated to annual GPP (Table 5.4). Figure 2 5.8 is a seasonal GPP plot from 2001 to 2010 with GPP in the unit of 0.0001 Kg C/(m *day). GPP change in this area from 2001 to 2010 has three stages: 2001 - 2002, 2003 - 2008, and 2009 - 2010. Highest annual GPP is found during 2001 to 2002. In 2003, a rapid drop of GPP occurred, very interestingly, corresponding to the starting time point of real estate development in Hangzhou. GPP remain stable during 2003 – 2008 except a relatively higher GPP value in 2007. Annual GPP dropped again in 2009. The overall trend is decreasing. From 2001 to 2010, mean GPP decreased by 23% while built-up area increased by over 80%. Figure 5.7 Land Cover Type in Hangzhou in 2001 59 Figure 5.8 Land Cover Type in Hangzhou in 2010 0.6 0.5 built-up 0.4 crop 0.3 tree 0.2 water 0.1 other 0 1 2 3 4 5 6 7 8 9 10 Figure 5.9 Land Cover Change in Hangzhou from 2001 to 2010 60 Table 5.4 Summary Statistics of Landsat-derived GPP in Hangzhou City Year 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Mean 5265.396 5116.157 4300.502 4427.953 4428.471 4434.384 4958.486 4354.339 3900.385 4049.857 Min 0 0 0 0 0 0 0 0 0 0 Max 13958.162 13682.897 10839.121 11417.204 11509.028 11811.452 13022.285 11554.476 10837.497 12100.703 Std. Dev. 3047.003 3156.753 2627.860 2848.209 2810.447 2857.593 3216.418 2978.814 2829.646 3031.548 2 Mean, Min, Max: 0.0001 kg C/(m *year) The TMIN_scalar and VPD_scalar are not the main cause of GPP change across the ten-year period (Figure 5.11 - 5.15). Even through the variation range became larger for these two factors, the mean scalar of each 16-day period of each year didn’t vary greatly. The climate change effects were not the primary cause of GPP decrease in the Hangzhou City from 2001 to 2010. The main change of GPP was because of urbanization process changed large area of cropland into built-up area. In 2001, cropland occupied more than 50% of the city area and build-up area was less than 20%. In 2010, the area cropland and built-up land was almost equal. The mean GPP decreased by 23% from 2001 to 2010 most likely due to urbanization process. Hangzhou City has the most expensive average house price ($3800 per square meter) in China. With such economic stimulus, further development of built-up area can be expected, together with decrease in GPP. In the analysis, daily minimum, maximum, average and dew point temperature are from ground observation data. The data collection is often based on repeated measurement at multiple time-points. It makes sure the high quality data are acquired. On the other hand, remotely sensed data have higher uncertainties like cloud contamination. Studies that are more quantitative in nature may need to perform cloud screening processes to reduce this influence. 61 Mean Daily GPP in Hangzhou City from 2001 to 2010 GPP (0.0001 Kg C/(m2*Day) 30 25 20 15 Landsat 10 5 0 1 353 702 1051 1400 1750 2099 2448 2797 Figure 5.10 GPP in Hangzhou City from 2001 to 2010 62 3147 3496 Daily Mean Temperature from 2001 to 2010 Temperature ℃ 40 30 20 10 0 -10 0 365 730 1095 1460 1825 2190 2555 Day Figure 5.11 Daily Mean Temperature from 2001 to 2010 Table 5.5 Summary of Temperature Change from 2001 to 2010 Year 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Mean 17.425 17.669 17.503 17.720 17.444 18.171 18.361 17.490 17.802 17.410 63 Std. Dev. 8.531 8.111 9.409 8.700 9.841 8.972 8.769 9.432 9.330 9.072 2920 3285 3650 Daily T_scalar for Tree 1.2 1 T_scalar 0.8 0.6 0.4 0.2 0 0 365 730 1095 1460 1825 2190 2555 2920 Day Figure 5.12 Daily T_scalar for Forest Land Cover from 2001 to 2010 64 3285 3650 Daily T_scalar for Crop 1.2 1 T_scalar 0.8 0.6 0.4 0.2 0 0 365 730 1095 1460 1825 2190 2555 2920 Day Figure 5.13 Daily T_scalar for Crop Land Cover from 2001 to 2010 65 3285 3650 Daily VPD_scalar for Tree 1.2 1 VPD_scalar 0.8 0.6 0.4 0.2 0 0 365 730 1095 1460 1825 2190 2555 2920 Day Figure 5.14 Daily VPD_scalar for Forest Land Cover from 2001 to 2010 66 3285 3650 Daily VPD_scalar for Crop 1.2 1 VPD_scalar 0.8 0.6 0.4 0.2 0 0 500 1000 1500 2000 2500 3000 Day Figure 5.15 Daily VPD_scalar for Crop Land Cover from 2001 to 2010 67 3500 4000 5.2.2 Comparison of Landsat GPP and MODIS GPP This subsection will compare Landsat-derived GPP with MODIS GPP. Figure 5.14 shows the GPP derived from two algorithms from 2003 to 2010 (23 periods per year). Mean daily GPP value was compared only in the area where MODIS GPP is not using fill value for built-up. In other words, the comparison is mainly on cropland and forest. Even though MOD 17 product also shows a decreasing trend in GPP, the MOD 17 calculated GPP has a higher value than landsat-derived GPP (Table 5.5). Figure 5.15 plots two algorithm’s result together, showing the relationship of Landsat GPP and MOD 17 GPP. The overall trend of two algorithms seems to fit. But comparing to Landsat-derived value, MODIS GPP has much higher values during growing seasons. In some years, EVI-derived GPP is nearly the half of MODIS GPP during the centermost periods of each year. MODIS GPP has more spikes pointing upward or downward, comparing to Landsat-GPP. The GPP algorithm in MOD 17 and LUE used in this study is based on same theory, calculated as: GPP =εmax * TMIN_scalar * VPD_scalar * FPAR * SWRad * 0.45 The main difference in the equation is FPAR. MOD 17 uses MOD 15 FPAR product as FPAR inputs, while this study adopts EVI as equivalent of FPAR. Wu et al. (2008) also found smaller EVI than NDVI in their modeling of GPP. They concluded that FPAR estimation based on NDVI-FPAR function is likely to represent the FPAR absorbed by all green vegetation, while EVI-FPAR function only represent the FPAR absorbed by leaf chlorophyll and EVI-GPP function is more linear related (Wu et al, 2008). As EVI is harder to get saturated in densely vegetated area, it is not as noisy as NDVI. In theory, whereas the NDVI is chlorophyll sensitive, the EVI is more responsive to canopy structural variations. The FPAR-EVI function has been examined in temperate forest and temperate grassland. 68 The results of those studies generally agree with the observation based on flux tower measurement. The vegetation in Hangzhou is consisting of forest and large area of cropland. However, two land cover types both present a lower GPP than MODIS GPP. Many studies found that MODIS GPP underestimate or overestimate GPP in many cases when comparing MODIS GPP to field flux tower observed GPP (Sjöström et al, 2011; Stagakis et al, 2007). In this case, the use of EVI as an equivalent of FPAR may underestimate the GPP because the vegetation structure in this area is not very densely vegetated so that using FPAR-EVI function didn’t respond well to chlorophyll, which is closely related to photosynthesis. Given the linear relationship of EVI and FPAR that has been reported by many studies, the coefficient “a” may need further analysis, for this specific study site. Many studies reported EVI-GPP relationship with great variation. The contamination of remotely sensed VI and different seasonality of EVI and GPP are two main issues: 16-day composite may still suffer from cloud contamination when there is no cloud-free acquisition during a 16-day period; during leaf expansion period, leaf photosynthetic capacity and leaf area change may not be coincident (Nagai et al., 2010). Studies should be carried out to observe GPP in the field with flux tower and compare the result with two algorithms addressed in this study: the difference of remotely sensed VI and ground-observed VI and the difference of remotely sensed GPP and ground-observed GPP need to be examined for the specific site. In the algorithm, TMIN_scalar and VPD_scalar are two factors having direct influence on the ε value. From 2001 to 2010, the daily mean temperature went up slightly (Figure 5.11) and seasonal temperature variation is greater during 2006 to 2007. During spring and winter, main limitation for ε is temperature. In the summer, vapor pressure deficit is the main limit (Figure 5.14-5.15). The MODIS Algorithm uses daily DAO with spatial resolution of 1.00°×1.25° as meteorological inputs, while this study used ground-station data. As GPP value is sensitive to meteorology inputs, higher temperatures may result in high vapor 69 pressure deficits and lower GPP. It is likely that the approach of adopting ground-station observation data may introduce uncertainty of GPP estimation to some extent in the urban-nonurban transition. However, as efforts were made to assess the possibility of deriving scale factors from remotely sensed data, the MODIS daily temperature product and evaporation product usually have a substantial number of missing values over a very large area. The same situation happens with 8-day composite products. MODIS monthly average products are less affected by missing data issue. Balance between spatial and temporal resolution is of primary concern when making decisions about which product to use. 70 Mean Daily GPP Derived from Two Algorithms from 2003 to 2010 50 GPP (0.0001 Kg C/(m2*Day) 45 40 35 30 25 MOD 17 20 Landsat 15 10 5 0 1 23 45 67 89 111 133 155 Figure 5.16 Mean Daily GPP Derived from Two Algorithms in Hangzhou City 71 177 30 25 Landsat 20 15 10 5 y = 0.527x + 3.3868 R² = 0.6211 0 0 10 20 30 40 50 MOD 17 Figure 5.17 The relationship of Landsat GPP and MOD 17 GPP 5.3 Method Derivation for Estimating EVI 5.3.1 Vegetation Phenology in Hangzhou City The typical time-series of built-up, tree and crop land cover types are shown in Figure 5.18-5.20. Built-up area shows seasonality because large ratio of tree and grass cover (Figure 5.18) with a lower EVI than the other two land cover types (Figure 5.19 – 20). Forest in the area has one season per year while the situation for crop is more complex (one season or two seasons). The human impact on cropland is significant, thus forecast or backcast on EVI of crops is unreliable. The uncertainty may come from the change of crop type or change of farming habits. In this study, only tree land cover is used for EVI estimation method development. 72 Figure 5.18 Time-series EVI for Built-up Area 73 Figure 5.19 Time-series EVI for Forest Land Cover Type 74 Figure 5.20 Time-series for Crop Land Cover Type 75 5.3.2 Function for Estimating EVI The F(t) function representing the basic trend of vegetation seasonality. Using the F(t) function to fit to time-series data for data from three windows, the R-square generally ranges from 0.4 to 0.77: within 120 observations, 15% are above 70%, more than 51% range from 0.6 to 0.7, the rest are below 0.6. Stepwise analysis was performed to construct EVI function with climate indicators (temperature and precipitation). Temperature and precipitation data were aggregate into 16-day average and 16-day summation, accordingly. The regression results indicate that the basic seasonal variation has been taken out by F(t) function. Using temperature and precipitation as independent variables to explain the interannual difference does not work well. One the other hand, the seasonal cycle has very high correlation with temperature (generally larger than 0.9). The main seasonal variation is attributed to temperature or radiation. Figure 5.21 to Figure 5.23 present the relationship of EVI and temperature (Ti, where i represents how many periods before EVI acquisition) using one cell as example. Figure 5.24 presents the relationship of EVI and precipitation (Pi, where i presents how many periods before EVI acquisition). The fitted-EVI reduced the noise of cloud contamination greatly and has a better correlation with temperature (Figure 5.27). As F(t) function or EVI-T regression generally underestimate EVI around peak value, temperature is the main factor affects vegetation growth, and T0 have highest correlation with EVI, equation (3) was established: 2 EVI_est = F1 * T0 + F2 * T0 + F3 where F1, F2, and F3 are one set of coefficients for one pixel, T0 is the 16-day period at the same 16-day of EVI acquisition. The result is shown in Figure 5.26. 76 0.6 0.5 0.5 MODIS EVI 0.7 0.6 MODIS EVI 0.7 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 10 20 T 0 (℃) 30 40 0 30 40 0.5 MODIS EVI 0.6 0.5 20 T 1 (℃) 0.7 0.6 MODIS EVI 0.7 10 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 10 20 T 2 (℃) 30 0 40 10 20 T 3 (℃) Figure 5.21 The Relationship of EVI and Temperature (T0 to T3) 77 30 40 0.6 0.5 0.5 MODIS EVI 0.7 0.6 MODIS EVI 0.7 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 10 20 T 4 (℃) 30 40 0 30 40 30 40 0.5 MODIS EVI 0.6 0.5 20 T 5 (℃) 0.7 0.6 MODIS EVI 0.7 10 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 10 20 T 6 (℃) 30 40 0 10 20 T 7 (℃) Figure 5.22 The Relationship of EVI and Temperature (T4 to T7) 78 MODIS EVI MODIS EVI 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 15 0 5 P1 15 10 15 P2 0.7 0.6 0.6 0.5 0.5 MODIS EVI 0.7 MODIS EVI 10 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 5 10 0 15 5 P3 P4 Figure 5.23 The Relationship of EVI and Precipitation (P1 to P4) 79 0.6 0.5 0.5 MODIS EVI 0.7 0.6 MODIS EVI 0.7 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 5 P5 10 0 15 0 5 10 15 10 15 0.5 MODIS EVI 0.6 0.5 P6 0.7 0.6 MODIS EVI 0.7 5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 5 P7 10 15 P8 Figure 5.24 The Relationship of MODIS EVI and Precipitation (P5 to P8) 80 0.6 0.5 0.5 EVI_fitted 0.7 0.6 EVI_fitted 0.7 0.4 0.3 0.2 0.4 0.3 0.2 0.1 0.1 0 0 0 10 20 30 40 0 10 T1 20 30 T2 Figure 5.25 The Relationship of Function-fitted EVI and Temperature (T1 to T2) 0.7 0.6 0.5 evi_fitted 0.4 evi_est 0.3 0.2 1 23 45 67 89 111 133 155 177 Figure 5.26 Estimated EVI from 2001 to 2010 81 199 221 40 Temperature at the same 16-day period of EVI acquisition has the best linear relationship with EVI. With the increase of time-period, the results shift from positively correlated to less correlated, then to negatively correlated. The EVI-precipitation relationship, on the other hand, shows a very low correlation. It indicates that precipitation is not the main limiting factor causing seasonal variation of tree growth in Hangzhou as Hangzhou has a humid subtropical climate with abundant rainfall all the year around. Other factors including soil moisture, relative humidity should be introduced into the estimation function to test their correlation to EVI. Figure 5.26 shows the estimated EVI and function-fitted EVI. The two EVI generally have high correlation (larger than 0.8). However, most difference occurs the peak valley time. The temperature factor along is not good enough to be used in the forecast or backcast. The quality of EVI may be another problem: data from the satellite suffers from cloud contamination even in a 16-day composite; the fitted function is limited by original data, which introduce further uncertainty in the later analysis. Vegetation grows in a suitable temperature. The lowest and highest temperatures in the Hangzhou area are not severe enough to stop tree growth for a long period. So during the aggregate of temperature, lower temperature was not ruled out. Generally, the estimated EVI was not accurate enough for GPP analysis due to the shifted peak and valley. T0 has better correlation than other periods, which is probably because T0 has better linear relationship with the effect of “growing degree day”. Growing degree day quantifies the cumulative effect of temperature. Further analysis could be carried out to examine phenology change based on growing degree day. Being a factor controlling vegetation stomata status, vapor pressure deficit was also examined with regression. Combination of VPD, temperature and precipitation all yields lower correlation than F(t) function. 82 Further analysis of adopting Gaussian function or other functions to derive EVI estimation functions is not conducted in this study. The original intention of deriving EVI estimation function is to using the curve and a single time-point EVI to derive EVI for the whole year. The ideal estimation is based on two parts: basic seasonal cycle of vegetation phenology and interannual variation that can be linked to climate indicators or other factors. To satisfy this requirement, factors play as the driving factors of interannual variability other than temperature should be identified or better estimation function based on temperature should be derived first, in order to yield more accurate estimation. In this case, if human disturbance is involved to a large extent, the activities are hard to quantify and thus not suitable for prediction. One finding during this part of study is that for many studies use VI as inputs, the cloud contamination issue can be reduced with suitable filtering method, if the studies are not carried out in environment highly disturbed by human activities. 83 6 Conclusion In this study, efforts are made to assess urbanization and its carbon consequences in the Yangtze River Delta. First, an empirical model was established based on MODIS GPP, population, GDP and built-up area data with two equations. The empirical work attempts to answer what is the impact of urban expansion on GPP through a balanced panel of 16×6=96 observations (one for each city-year from 2003 to 2008 within the 16 cities located in Yangtze River Delta). Second, GPP in the Hangzhou city area was calculated based on the light-use efficiency method with MODIS EVI, Landsat image, NOAA ground-station climate data, and GLDA radiation data. The land use and land cover change was assessed from 2001 to 2010 in this area. The Landsat-derived GPP was compared to MODIS GPP. Thirdly, the method of estimating EVI based on vegetation phenology and climate indicators was assessed in order to estimate GPP in a data-sparse environment when VI products with high temporal resolution were not available (in 1970s, only one to two scenes of Landsat image can be used to derive vegetation indices every year). According to the empirical model estimations in the Yangtze River Delta, urban expansion is significantly correlated to decrease of both the annual GPP sum and GPP density. Column 1 in Table 5.2 shows that 1) a one-percent increase in built-up area (i.e. the built-up area increased from 100 km2 to 101 km2) is significantly correlated to a 1.066-percent decrease in annual GPP sum; a one-percent increase in population is significantly correlated to the by a 1.176-percent decrease of annual GPP sum; Column 2 in Table 5.2 shows that a one-percentage increase in built-up area ratio (i.e. the land area for a city is 1000 km2, within which the built-up area ratio increased from 10 percent (100 km2) to 11 percent(110 km2) ) is strongly correlated to a 15.4-percent decrease in annual GPP sum. The GPP density equation yields similar results. Another interesting finding is that capital cities has more GPP sum and 84 density compared with non-capital cities. The possible reason is that those cities have more funding and focus more on environment such as a more strict regulation for pollute-intensive industry formation. Capital cities are generally key tourism regions: they are more concerned about increasing green space and attracting tourists. Finally, the estimation of year fixed effects predicts a decrease trend of GPP sum and density across time. In the Yangtze River Delta, the efficiency of land use is one major issue. In 1992, the development strategy led to rapid establishment of new development zones. In 1991, there were only 17 developing zones in the area. By 1993, this number increased to 74. During this process, large areas of agricultural land were enclosed blindly and standing idle after enclosure. In some cities, the actual land enclosed is three times higher than approval (Liu, 2005; Yang, 2001). Industrial uses of land were not wisely planned in many cases. The first part of this study shows the decreasing GPP with further urbanization in the Yangtze River Delta quantitatively. If this process continues to go on, primary production will continue to decline. This may increase the possibility of land degradation with heavier use of fertilizer. To achieve a sustainable future, land use in the Yangtze River Delta should be planned more wisely, and the structure of development should also be changed. To get rid of the dilemma of development and environment protection, an intensively sustainable development policy must be carried out to support a smart growth of cities in the area. In the Hangzhou city area, land cover type has changed dramatically since 2001. In 2001, cropland occupied more than 50% of the city area and build-up area was less than 20%. In 2010, the areas of cropland and built-up land were almost equal. The mean GPP decreased by 23% from 2001 to 2010 along with the urbanization process. It is very likely that the GPP will continue decreasing in the future with further urban expansion. Comparing to MODIS GPP, Landsat GPP has an overall similar trend. But MODIS GPP has much higher values during growing seasons. As many studies have reported that MODIS GPP tend to 85 underestimate GPP comparing to ground observation GPP, the landsat-derived GPP is likely underestimated. The main difference in the algorithm is that Landsat-derived GPP adopts FPAR-EVI function. Landsat-derived GPP may underestimate GPP as EVI is more responsive to canopy structural variations and does not respond well to chlorophyll. In order to make further conclusion, studies should be carried out to observe GPP in the field. During the simulation of EVI, regression results show that using F(t) function as seasonal cycle and try to explain interannaul difference with temperature and precipitation does not work well for the study sites. Temperature at the same 16-day period of EVI acquisition has the best linear relationship with EVI (Figure 5.21 - 5.23) while the EVI-precipitation relationship shows a very low correlation. The estimated EVI’s shifts in peaks and valleys comparing to function-fitted EVI suggest that using estimated EVI is not good enough for GPP analysis. To satisfy the need to estimate VI in a data-sparse environment, factors play as the drivers of interannual variability other than temperature should be identified first or better simulation function should be derived. If human disturbance is the main reason for interannual difference, this method is not suitable for backcast, as human activities are hard to quantify. On the other hand, the study finds filtering process over noisy VI outputs may reduce cloud contamination. However, caution is needed during filtering in case the real change is eliminated. This study demonstrates that urbanization processes in the Yangtze River Delta has a negative correlation with gross primary production. Landsat imagery is suitable for urbanization analysis. During the estimation of GPP, FPAR-EVI function needs to be examined with ground-truth. Derivation of other functions to better simulate interannual variation of VI will contribute to better estimation of VI. However, the factors could involve human disturbance, which is hard to quantify. Additionally, satellite products themselves contain inevitable errors, which also add difficulties to analysis. 86 Further studies could investigate the demand change of primary production in the Yangtze River Delta to better understand the balance and budget of GPP. Future studies could also investigate how to derive VPD based on satellite data with better spatial and temporal resolution. More ground-measured FPAR-EVI relationship based on different vegetation types and geographic areas would certainly provide insights on GPP analysis. Grassland VI prediction with temperature and precipitation yields results agree well with actual value (Iwasaki, 2009). The prediction resolution is relatively coarse (i.e. 0.25°×0.25°). Vegetation phenology is one important tool to improve VI forecast or backcast with higher spatial or temporal resolution. For satellite products with substantial errors due to cloud contamination, suitable filtering method would reduce the cloud effects. In future similar projects, filtered vegetation indices should be adopted. Also, nighttime light generated remote sensing products for GDP, population distribution can be used to assess spatial distribution of GPP and those factors. 87 Appendices 88 ERDAS Script for Changing Fill Values in Raster File to Zero COMMENT "Generated from graphical model: Change Fill Value"; # # set cell size for the model # SET CELLSIZE MIN; # # set window for the model # SET WINDOW UNION; # # set area of interest for the model # SET AOI NONE; # # declarations # Integer RASTER n1_2002_185_clip FILE OLD NEAREST NEIGHBOR AOI NONE arg1; Integer RASTER n3_jiaxing FILE DELETE_IF_EXISTING USEALL ATHEMATIC 16 BIT SIGNED INTEGER arg2; # # function definitions # n3_jiaxing = CONDITIONAL { ($n1_2002_185_clip >= 30000) 0 , ($n1_2002_185_clip < 89 30000) $n1_2002_185_clip } ; QUIT; 90 IDL Script for Coverting GeoTiff File to Flat Binary File pro tif2bin files = dialog_pickfile(filter='*.tif',/read,/multiple,title='Please input images') num=n_elements(files) path = 'J:\EVI \' FOR i =0, num -1 DO BEGIN tmp = read_tiff(files[i],geotiff = tifinfo) tmp2 = float(tmp)/10000.0 newfile = path + file_basename(files[i],'tif') OPENW, 1, newfile ; Write the data into the file: ; FIX() function is used to change the file format to 16bit signed integer!! WRITEU, 1, tmp2 ; Close file unit 1: CLOSE, 1 ENDFOR print, 'Congratulations!!!' end 91 Bibliography 92 Antrop, M.. 2004. Landscape change and the urbanization process in Europe. Landscape and Urban Planning. 67: 9-26. Ardö, J., M. Mölder, B. El-Tahir and H. Elkhidir. 2008. Seasonal variation of carbon fluxes in a sparse savanna in semi arid Sudan. Carbon Balance and Management. 3: 7 Bounoua, L., Safia, A., Masek, J., PetersLidars, C., and Imhoff, M. L. (2009). Impact of urban growth on surface climate: A case study in Oran, Algeria. J. Appl. Meteor. Clim., 48, 217-231. Cervero, R., and Day, J. 2008. Suburbanization and transit-oriented development in China. Transport Policy, 15(5): 315-323. Champion, T., 2001. Urbanization, suburbanisation, counterurbanisation and reurbanisation. In: Paddison, R. (Ed.), Handbook of Urban Studies. Sage, London, pp. 143 – 161. Chen, S., Zeng, S., and Xie, C.. 2000. Remote Sensing and GIS for Urban Growth Analysis in China. Photogrammetric Engineering and Remote Sensing. 66 (5): 593 – 598. Curtis, O., Hanson, P., Bolstad, P., Barford, C., Randolph, J., Schmid, H., Wilson, L.. 2002. Biometric and Eddy Covariance Based Estimates of Annual Carbon Storage in Five Eastern North American Deciduous Forests. Agricultural and Forest Meterology. 113: 3 – 19. Dell, A., and Gamba, P.. 2001. Detection of urban structures in SAR images by robust fuzzy clustering algorithms: The example of street tracking. IEEE Trans. Geosci. Remote Sensing 39 (10): 2287–2297. Eklundh, L., and Jönsson, P., 2010, TIMESAT 3.0 - Software Manual. Lund University, 74 pp. Eklundh, L. and Jönsson, P., 2009, Timesat 3.0 Software Manual, Lund University, Sweden. Jönsson, P. and Eklundh, L., 2002, Seasonality extraction and noise removal by function fitting to time-series of satellite sensor data, IEEE Transactions of Geoscience and Remote Sensing, 40, No 8, 1824 – 1832. Fensholt, R. I., Sandholt and M.S. Rasmussen. Evaluation of MODIS LAI, fAPAR and the relation between fAPAR and NDVI in a semi-arid environment using in situ measurements. 2004. Remote Sensing of Environment 91: 490–507 Fu, S., and Gao, G.. 1989. Jing-Jin-Tang Atlas of Ecological Environment. Science Press. Beijing. 316p. G. Grigera, M. Oesterheld. Forage Production monitoring system : parameterization of fPAR and RUE. 2006. Poster presented at "Global Vegetation Workshop". University of Montana. IFEVA, Facultad de Agronomia, Universidad de Buenos Aires. http://www.ntsg.umt.edu/VEGMTG/posters/Griger_VegWorkshop_Missoula2006.ppt Gallo, K. P., McNab, A. L., Karl, T. R., Brown, J. F., Hood, J. J., and Tarpley, J. D.. 1993. The use of a vegetation index for assessment of the urban heat island effect. International Journal of Remote Sensing. 14: 2223- 2230. Göckede, M., T. Foken, M. Aubinet, M. Aurela, J. Banza and C. Bernhofer et al.. 2008. 93 Quality control of CarboEurope flux data—Part 1: coupling footprint analyses with flux data quality assessment to evaluate sites in forest ecosystems, Biogeosciences 5: 433–450 Goward, S., Tucker, C., Dye, D.. 1985. North American Vegetation Patterns Observed with the NOAA-7 AVHRR. Vegetatio. 64: 3 – 14. Güneralp, B., and Seto, K.. 2008. Environmental Impacts of Urban Growth from an Integrated Dynamic Perspective: A Case Study of Shenzhen, South China. Global Environmental Change. 18: 720 – 735. Haack, B., Bryant, N., and Adams, S., 1987. Assessment of Landsat MSS and TM Data for Urban and Near-Urban Landcover Digital Classi. cation. Remote Sensing of Environment. 21: 201–213. He, J., and Zhuang, D..2006. Spatial and Temporal Land Use Pattern and Environment Impacts in Yangtze River Delta. Geographical Research. 25 (13): 388 – 396. (in Chinese) Imhoff, M., Bounoua, L., Ruth, D., Lawrence, W., Stutzer, D., Tucker, c., and Ricketts, T.. 2004. The Consequences of Urban Land Transformation on Net Primary Productivity in the United States. Remote Sensing of Environment. 89: 434 – 443. Jensen, J., and Cowen, D.. 1999. Remote Sensing of Urban/Suburban Infrastructure and Socio-Economic Attributes. Photogrammetric Engineering and Remote Sensing. 65 (5): 611-622. Ji, M., and Jensen, J., 1999. Effectiveness of subpixel analysis in detecting and quantifying urban imperviousness from Landsat Thematic Mapper imagery, Geocarto International. 14(4): 31–39. Jing, X. 2007. Suburbanization Trend of Large Cities in the Yangtze River Delta (in Chinese). Zhejiang Economy. 11: 17 – 19. Jönsson, P. and Eklundh, L., 2004, Timesat - a program for analyzing time-series of satellite sensor data, Computers and Geosciences, 30, 833 – 845. Kanniah, K.D., J. Beringer, L.B. Hutley, N.J. Tapper and X. Zhu. 2009. Evaluation of collections 4 and 5 of the MODIS Gross Primary Productivity product and algorithm improvement at a tropical savanna site in northern Australia, Remote Sensing of Environment 113: 1808–1822 Lieth H., Whittaker R.. 1975. Primary Productivity of the Biosphere. New York. Springer-Verlag. Lindroth, A., F. Lagergren, M. Aurela, B. Bjarnadottir, T. Christensen and E. Dellwik et al.. 2008. Leaf area index is the principal scaling parameter for both gross photosynthesis and ecosystem respiration of Northern deciduous and coniferous forests, Tellus Series B-Chemical and Physical Meteorology 60: 129–142 Liu, h. 2005. Situation of land resources in yangtze river delta and counter measures for sustainable utilization. Chinese journal of agricultural resources and regional planning. 26(4): 9 -13. 94 Lu, D., Xu, X., Tian, H., Moran, E., Zhao, M., Running, S.. 2010. The Effects of Urbanization on Net Primary Productivity in Southestern China. Environmental Management. 46: 404 – 410. Ma, L. 2004. Economic reforms, urban spatial restructuring, and planning in China. Progress in Planning. 61: 237 – 260. Monteith, J.. 1972. Solar Radiation and Productivity in Tropical Ecosystems. Joural of Applied Ecology. 9: 747 – 766. Monteith, J.. 1977. Climate and Efficiency of Crop Production in Britain. Philosophical Transactions of Royal Society of London, Ser. B: 277 – 294. Nagai, S, Saigusa, N, Muraoka, H, Nasahara, KN. 2010. What Makes The Satellite-Based Evi-Gpp Relationship Unclear In A Deciduous Broad-Leaved Forest? Ecological Research, 25(2), 359-365. Ni, J.. 2004. Estimating Net Primary Productivity of Grasslands from Field Biomass Measurements in Temperate Northern China. Plant Ecology. 174: 217 - 234. Niu, Z., Wang, C.. Basis and Application of Carbon Cycle Remote Sensing. Science Press. Beijing. 2008. 286p. Pan, X. and Zhao, Q.. 2007. Measurement of Urbanization Process and the Paddy Soil Loss in Yixing City, China between 1949 and 2000. Catena. 69: 65-73. Pant, M. 2008. Response of Vegetation Phenology to Rainfall timing in the Sahel 1982 to 2004. Thesis Submitted to the International Institute for Geo-information Science and Earth Observation. Rottensteiner, F., Briese, Ch., 2002. A new method for building extraction in urban areas from high resolution LIDAR data. ISPRS. Photogrammet. Comput. Vision, Graz, Austria,9–13 September, pp. A-295 ff. Running, S., Nemani, R., Heinsch, F., Zhao, M., Reeves, M., and Hashimoto, H.. 2004. A Continuous Satellite-Derived Measure of Global Terrestial Primary Production. BioScience. 54 (6): 547 – 559. Running, S., Nermani, R.. 1988. Relating Seasonal Patterns of the AVHRR Vegetation Index to simulated Photosynthesis and Transpiration of Forests in Different Climates. Remote Sensing of Environment. 24: 347 -367. Running, S., R. R. Nemani, F. A. Heinsch, M. Zhao, M. Reeves, and H. Hashimoto. 2004. A continuous satellite-derived measure of global terrestrial primary production. BioScience. 54(6): 547-560. Seto, K., and Fragkias, M.. 2005. Quantifying Spatiotemporal Patterns of Urban Land-use Change in Four Cities of China with Time Series Landscape Metrics. Land Ecology. 20: 871 – 888. Sicular, T. 1985. Chian’s grain and meat economy recent development and implications for trade. American Journal of Agricultural Economics. 67: 1055 – 1062. 95 Song, C. 2005. Spectral mixture analysis for subpixel vegetation fractions in the urban environment: How to incorporate endmember variability? Remote Sensing of Environment, 95, 248?263. Stavros, s., Markos, N., Levizou, E., and Kyparissis, A.. 2007. Envisat Symposium 2007. Montreux, Switzerland. Streutker, D.. 2003. Satellite-measured Growth of the Urban Heat Island of Houston, Texas. Remote Sensing of Environment. 85: 282 – 289. Tans, P., Fung, I., and Takahashi T.. 1990. Observational Constraints on the Global Atmospheric CO2 Budget. Science. 247: 1431 – 1438. Tetens, O., 1930. Uber einige meteorologische Begriffe. z. Geophys. 6:297-309. Tucker, C., Townshend, J., and Goff, T.. 1985. African Landcover Classification Using Satellite Data. Science. 227: 369 – 374. Turner, D., Urbanski, S., Wofsy, S., Bremer, D., Gower, S., &and Gregory, M.. 2003. A cross-biome comparison of light use efficiency for gross primary production. Global Change Biology. 9: 383–395. Wang, Z., Xiao, X., Yan, X.. 2010. Modeling Gross Primary Production of Maize Cropland and Degraded Grassland in Northestern China. Agricultural and Forest Meteology. 150: 1160 – 1167. Welch, R.. 1980. Monitoring Urban Population and Energe Utilization Patterns from Satellite Data. Remote Sensing of Environment. 9: 1 - 9. Weng, Q.. 2002. Land use change analysis in the Zhujiang Delta of China using satellite remote sensing, GIS and stochastic modelling. Journal of Environmental Management. 64 (3): 273-284. Wu, W., Wang, S., Xiao, X., Yu, G., Fu, Y., and Hao, Y.. 2008. Modeling gross primary production of a temperate grassland ecosystem in Inner Mongolia, China, using MODIS imagery and climate data. Sci China Ser D-Earth Sci. 51 (10): 1501 – 1512. Xiao, J., Shen, Y., Ge, J, Tateishi, R., Tang, C., Liang, Y., and Huang, Z.. 2006. Evaluating Urban Expansion and Land Use Change in Shijiazhuang, China, by Using GIS and Remote Sensing. Landscape and Urban Planning. 75: 69 – 80. Xiao, X., Zhang, Q., Saleska, S., Hutyra, L., De Camargo, P.,Wofsy, S., et al. 2005. Satellite-based Modeling of Gross Primary Production in a Seasonally Moist Tropical Evergreen Forest. Remote Sensing of Environment. 94: 105–122. Xu, C., Liu, M., An, S., Chen, J., and Yan, P.. 2007. Assessing the Impact of Urbanization on Regional Net Primary Productivity in Jiangyin County, China. Journal of Environmental Management. 85: 597 – 606. Xu, H., Wang, X., and Xiao, G.. 2000. A Remote Sensing and GIS Integrated Study on Urbanization with Its Impact on Arable Lands: Fuqing City, Fujian Province, China. Land Degradation and Development. 11: 301 – 314. 96 Yang, G. 2001. The process and driving f orces of change in arable2land area in The yangtze river delta during the past 50 years. Journal of Natural Resources. 16 (2): 121 127. Yu, D., Shao, H., Shi, P., Zhu, W., and Pan, Y.. 2009. How Does the Conversion of Land Cover to Urban Use Affect Net Primary Productivity? A Case Study in Shenzhen City, China. Agricultural and Forest Meteorology. 149: 2054 – 2060. Zhao, T. 2007. Changing Primary Production and Biomass in Heterogeneous Landscapes: Estimation and Uncertainty Based on Multi-Scale Remote Sensing and GIS Data. Dissertation submitted to the University of Michigan. Zhao, T., Brown, D., Bergen, K.. 2007. Increasing Gross Primary Production (GPP) in the Urbanizing Landscapes of Southeastern Michigan. Photogrammetric Engineering and Remote Sensing. 73 (10): 1 – 9. 97