Mimika Land Rights Project - West Papua
Prepared by Mike Alonzo (GISCorps Volunteer) and Nabil Ahmed (Forensic Architecture)
West Papua is the remote Melanesian province of Indonesia that has an abundance of gas, mineral and timber resources (Fig 1). Together with Papua New Guinea, in terms of biodiversity, West Papua contains 5% of the world’s tropical high-biodiversity terrestrial ecosystem. At the same time West Papua is a conflict zone with heavy Indonesian military presence, as there is a decade’s long independence movement waged by Papuans. The military has a long history of dealing heavy-handed with political protest and civilians alike and there are numerous cases of well-documented human rights violations that amounts to a slow genocide.
Figure 1: The study area in Western Papua
West Papua remains a region under close supervision by the military and a difficult area to operate in for human rights workers, environmental activists and journalists. The Grasberg mine operated by Freeport-Rio Tinto is located in Amungme territory in the Carstenz range. It is considered to hold the planet’s largest combined reserve of copper and gold. Copper from the Grasberg mine is the cheapest in the world.
The Grasberg mine has been dumping tailings into the Ajkwa River deposition area downstream in Komoro territory since operation in the 1980’s and is the site of an environmental disaster. At the same time, Freeport-Rio Tinto have a documented history of financing Indonesian military presence in West Papua, which activists claim is supporting a slow genocide.
Nabil Ahmed (Forensic Architecture and Earth Sensing Association) has launched an investigation into documenting the Grasberg mine environmental disaster on behalf of Octavianus Mote, Papuan negotiator and in collaboration with TAPOL, the human rights campaign organization for Indonesia based in the UK. The spatial analysis conducted by GISCorps Remote Sensing Specialist is a crucial preliminary element of the investigation, especially for an area that is extremely difficult to reach as a zone of conflict.
Objectives and Deliverables (i.e., tasks for the GISCorps Remote Sensing Specialist)
Objective #1: Map vegetation disturbance in Western Papua that can be attributed to either acid mine drainage from the Freeport-Rio Tinto mine or to population pressures (clearing for agriculture, urbanization) associated with increased mining activity.
- Deliverable 1: Animated visualization of land cover change from 1987 – 2014.
- Deliverable 2: Dated disturbance map from 1987 – 2014.
- Deliverable 3: Landscape complexity map illustrating the number of “states” of every pixel from 1987 – 2014.
Objective #2: Quantify the growth of the acid mine drainage area (ADA) by year and correlate with existing mine production data.
- Deliverable 1: Amount of land disturbed by year
- Deliverable 2: Visual comparison of disturbance with mine tailing production data
All data processing was completed in Matlab. Code will be made available by Mike Alonzo after publication of an associated journal article.
Overview: The basic plan was to create an algorithm which could track land cover change (anywhere in the world, potentially) based on change in any spectral index including but not limited to the Normalized Difference Vegetation Index (NDVI). The work described in this report was heavily inspired by the USDA Forest Service’s Landtrendr project. An adapted algorithm was required, however, because of the extreme cloudiness of the Western Papua region.
All cloudy areas were masked using the cfmask band (which only masks the most obvious clouds) and the 199 images were clipped to our study area extent. Significant clouds and haze remained. (Fig 2)
All reflectance images were transformed to NDVI which responds to the presence and density of vegetation. NDVI is calculated as (Near Infrared – Red) / (Near infrared + Red). We now have a data cube that is 2000 pixels (60 km north – south), 1300 pixels (39 km east – west) and 199 pixels in depth (1 for each date). (Fig 3)
If one clicks on a pixel that falls in the ADA, wanting to see the NDVI values over time, it will look something like this: (Fig 4)
If one clicks on stable forest, the NDVI values will look something like this: (Fig 5)
The goal now is to develop an algorithm that can fit several lines in a piecewise fashion to this rather noisy data in order to show: (1) the number of “states” of this pixel (e.g., 3 states could be stable forest à ADA disturbance event à stable ADA); (2) the type of disturbance (e.g., urbanization vs. ADA); and (3) the date of the disturbance. All of this must happen while avoiding “distraction” from the many data points that are bad because of cloud/haze intrusion.
The basis for the chosen algorithm is Matlab code called Shape Language Modeling (SLM) which was written by John D’Errico. This code makes use of non-linear optimization to simultaneously estimate the dates of a state change (e.g., change from stability to decline) as well as the associated NDVI value at that point. Figure 4 shows the result of this code visually using the red line for ADA disturbance. Figure 5 shows the single line for undisturbed forest. Figure 6 shows the process of urbanization. This is different from the ADA disturbance because it is more gradual (i.e., the slope of the red line is not as steep).
We ran this code for all 2.6 million pixels in the study area. This is very processing intensive and could literally take a month to run on a single processor. Luckily, many folks in the remote sensing community have access to some kind of multicore server. For this analysis, we used 32 processors in parallel and completed the study area in approximately 1.5 days.
Disturbance dates were assigned based on a change from higher than 0.40 NDVI to below 0.40 NDVI and a slope of the cross-over segment that was adequately high (a means for separating the somewhat ill-defined processes of “disturbance” and “decline”).
The year-by-year increase in the ADA was simply calculated as the number of pixels disturbed per year. This does include urbanized pixels as well but that number is fairly small and could be separated out simply if required. One pixel = 900 m2.
The land cover change animation: This deliverable was dealt with separately. Here we wanted to produce one relatively cloud-free composite image per year to add to an animation. To this end, we selected the pixel-date in a given year for each pixel location. “Best” was defined as smallest band1:band4 (blue to near infrared) ratio because high blue reflectance is broadly indicative of atmospheric contamination.
Figure 2: Study area during a particularly cloudy year. Solid green areas are the cfmask but note how many additional clouds remain unmasked
Figure 3: Left: NDVI image from 1987. Right: NDVI image from 2012. Red is high NDVI (~0.70 – 0.90) and blue is low NDVI (~0.10 – 0.30)
Some of the typical “pixel trajectories” (aka land cover histories) are shown in Figures 4,5, and 6.
Here is the map of estimated disturbance date. Black areas were not disturbed. (Fig 7)
Here is the map of landscape complexity. The interpretations of the complexities (with varying degrees of confidence are):
1 = stable forest (or less likely, stable river)
2 = decline/disturbance à stable ADA or urban area
3 = stable forest à disturbance à stable ADA
4 and 5: Probably complex, braided river dynamics possibly including disturbance and vegetation regrowth. (Fig 8)
Figure 4: NDVI change (translated to 0-100) over time in the ADA along with piecewise best-fit line (to be discussed later). This indicates a disturbance date of about 1998.
Figure 5: NDVI change for growing, undisturbed forest. The many instances of low NDVI can be considered noise related to atmospheric contamination (mostly clouds).
Figure 6: Urbanization characterized by a somewhat gradual loss of vegetation followed by a period of stability.
Figure 7: Disturbance date map
Figure 8: Landscape complexity. Or, how many things have happened to a given pixel. Most of the ADA has complexity of 3 which means it was first forested, then disturbed, and is now “river”.
Figure 9: Background (and left axis) is data from WALHI report on tailings production per year. Foreground (and right axis) is yearly growth of ADA. We hypothesize that after about 2000, the tailings were simply reaching the sea, thus the slowed growth of the ADA. This has not been rigorously confirmed.