2.1.

Materials

2.1.1.

Landsat

In this study, four scenes of multispectral (seven bands) Landsat-5 Thematic Mapper (TM) images have been used.²⁹ The images selected for this work correspond to the dates 13-4-2006, 02-7-2006, 19-08-2006, and 11-9-2006, and cover two areas of about $15\mathrm{km}\times 15\mathrm{km}$ (spatial resolution of 20 m) of the 197-031 and the 198-31 path-rows. They thus provide a set of images that is very representative of the main phenological cycle in this Mediterranean region. These two areas, in the Ripollès and Penedès counties, (Fig. 2, pink squares) are located in Catalonia, a region of approximately $32,000{\mathrm{km}}^2$ situated in the northeast Iberian Peninsula at the extreme southwest of Europe (map location Fig. 2). The main landscape of the Ripollès region is deciduous forest, while the Penedès region has a vineyard landscape. Landsat images were chosen because they are the remote sensing images most used for research purposes in different applications,³⁰ and they are also used in a large variety of applied projects.³¹

Fig. 2

Study regions: the Ripollès and Penedès areas for Landsat and MODIS are shown as pink rectangles and the Vallès and Selva areas for CASI images are shown as yellow rectangles.

2.2.

Image Processing

Different image processing methods were applied for each type of image depending essentially on the previous image processing level and on their particular characteristics.

Landsat images were acquired at the L1G processing level, and were then georeferenced and radiometrically corrected with MiraMon GIS software³⁴ according to the methodologies detailed in Refs. 35 and 36, and applied to the SatCat server.³⁷ The optical reflectance range of the corrected images extends from 0 to 10,000, representing the % of reflectance*100, and so numbers were written with two decimal places in a short integer. This factor and type of data (signed short integer) are also used to represent the possible range of ground temperatures in °C (derived from the radiometric correction) of the TM thermal band. All images use a $-999$ value to represent nodata regions, which in this case are normally caused by sensor problems or by post-processing artifacts. Defining the nodata value as $-999$ is not arbitrary because nodata values that are defined too far from the valid data range can negatively affect normalization compression procedures.³⁸ Consistent treatment of nodata values is one of the keys to performing correct image processing and subsequent data analysis.

However, CASI images were acquired geometrically georeferenced using the SISA system developed by ICC.³⁹ The flight orientation was SW-NE, which resulted in the geometrically corrected image being very large ($18,366\mathrm{columns}\times 11,918\mathrm{rows}$), and therefore the CASI images were rotated $-32.5\mathrm{degrees}$ with respect to the north projection in order to reduce the large amount of nodata pixels. The value range, data type and nodata definition in the CASI images were unified to match those of the Landsat images, also using the MiraMon software.

Finally, the MODIS images were downloaded from the Warehouse Inventory Search Tool⁴⁰ as a georeferenced reflectance product. It was only necessary to clip them to the study regions already detailed in the previous section and unify the nodata values.

2.3.

Lossy Compression

The reflectance product obtained was compressed at a wide range of different quality levels with the following compression ratios (CR): $2.5:1$, $5:1$, $10:1$, $50:1$, $100:1$, $200:1$, and $400:1$ (from soft to hard compression). Two methodologies, based on the JPEG2000 lossy compression procedures, were evaluated in this work:

Both compression procedures were applied using the Kakadu software,⁴¹ which is a complete implementation of the JPEG2000 standard, Part 1, completed with a great deal of Part 2 and 3.⁴² The BOÍ software,⁴³ an implementation of the JPEG 2000 (Part 1) standard, was used for calculating the PSNR quality compression indicator in reference to noncompressed images.

2.4.

Geostatistics

Geostatistics is a subset of statistics specialized in the analysis, interpretation, and inference of geographically referenced data.⁴⁴ It was initially developed by Matheron⁴⁵ and defined as a set of methods for studying the spatial distribution of a variable from a statistical approach in order to estimate the corresponding unknown value of a particular location or simulate its variability. When this variation has a spatial structure, the regionalized variable theory makes it possible to take spatial properties into account following a stochastic approach, and assumes a constant local mean and a stationary variance of the differences between regions separated by a given distance and direction.⁴⁶ The variogram (also referred to as semivariogram) is based on these previous assumptions and, as defined in Eq. (1), it plots the variance function depending on sample distance separation and, consequently, the spatial patterns analyzed. This expression implies, for remote sensing images, that pixel values variance is only a function of distance and direction between pixels; then, there is no dependence on their particularly locations and, in statistics sense, this variance pattern is stationary for all study regions. Therefore, the variogram describes the pattern distribution of image spatial resolution, and it can identify a study region, a subscene of remote sensing image in this work.

The sampled variogram can be modeled and fitted by a continuous function to identify structural parameters that characterize the spatial pattern for any quantitative variable distribution, including, of course, those of remotely sensed images, as carried out here in this study and in other previous studies⁴⁷ (Fig. 3). These parameters are:

Fig. 3

Variogram structure parameters.

For each specific geostatistical application, the key role parameter depends upon the focus of a particular spatial pattern analysis. For example, range is the most suitable for autocorrelation analysis, nugget for determining stochastic noise and sill for variability’s measurements.

The total sill was chosen as the main quality measure shown in Sec. 3 for analyzing alterations in the spatial patterns at studied compression ratios and with different methodologies. Different selected plots were used for three sensors for several dates, bands, and regions, which made it possible to determine the complete variogram structure by comparing the theoretical plot in Fig. 3 with the experimental plots in Sec. 3.

2.5.

Parallel Solution

Most variogram analyses are used to explore the spatial patterns of irregularly distributed samples.⁴⁸ Variogram analysis is usually a previous step to interpolation (kriging), which is a process with a high computational cost.⁴⁹

Nevertheless, the variogram analysis of remotely sensed images has higher computing demands to those of irregularly distributed samples. Indeed, a large number of pixels are involved even in geographically small scenes with a high or medium spatial resolution, and therefore a very large number of pairs are necessary to obtain the variogram. Table 1 compares four examples of the amount of data involved and the time necessary for executions running on a single, normal PC (Intel(R) Xeon (R) 3.0 GHz and 1 GB of RAM).

Table 1

Comparison of processing time of variogram analyses between two irregularly distributed point measurements and remote sensing images. The sparse irregular distribution corresponds to a network of weather stations, and the dense, irregular, distribution to a GPS network of altimetry measurements, both detailed in Ref. 44; the remotely sensed image corresponds to two (Landsat and CASI) of the three examples used in this work.

In this work, a large number of spatial pattern analyses have been executed. All spectral bands (the number of bands for Landsat, CASI, and MODIS are 7, 72, and 2, respectively) on several scene dates (4, 2, and 73) and in two different study regions were analyzed. In addition, the images were compressed at several compression ratios (eight in all cases) and compared with noncompressed images. This implied 504 analyses for Landsat, 2592 for CASI, and 2628 for MODIS. Performing all these analyses would be a very time-consuming process; a single CASI analysis takes 57,727 s (see Table 1). Therefore, a distributed environment was the most suitable solution for a complete and efficient study.

In order to reduce the execution time, the authors carried out a parallel implementation of the variogram analyses. The master/worker programming paradigm was used in which the master process schedules, distributes, and coordinates the tasks executed by the worker processes. The code was written in ANSI-C language, including MPI (message passing interface) as a message-passing library (version 1.2.7) that is used by the master process to communicate with the worker processes. Since MPI has become a standard de-facto message passing library,⁵⁰ this solution ($\mathrm{ANSI}\text{-}\mathrm{C}+\mathrm{MPI}$) guarantees portability to different computer platforms.

The distributed load design⁵¹ is not specifically defined and optimized for the characteristics of the images used in this work. However, it is a flexible design, tested with a wide range of image dimension that maintains a good balance and scalability suitable for other studies with a wide range of image dimensions and cluster properties. The scalability solution (adapting point samples analysis from Ref. 49 to remote sensing image analysis in the present work) can be achieved by defining a relatively small load unit: of two image rows of pixels. Each worker analyzes all possible combinations of nonrepeated pairs between two rows of the image and returns the partial variance and distance of the pairs involved to the master. The master distributes tasks (rows of data) to workers, accumulates partial results and manages input (remotely sensed image) and output (resulting variogram) tasks.

3.1.

Variogram Analysis

As explained above, lossy compression mechanisms modify the original data information. Figure 4 illustrates the effects of compression ratios on quality visualization of the images. In this figure, the central image is a noncompressed CASI example, the top image sequence corresponds to three BIFR compressed images; from left to right, it is shown increasing compression ratios and decreasing quality. The bottom image sequence corresponds to the same sequence improved using 3d-DWT techniques. This composition summarizes most of the results presented in this section, and the following tables and figures quantify these effects in a geostatistical sense using the variogram spatial pattern analysis.

Fig. 4

The effects of compression on a CASI image fragment from the Selva region: noncompressed in the middle, upper sequence of BIFR compression ratios: $2.5:1$, $10:1$, and $400:1$. The images in the bottom sequence were compressed at the same ratios but the 3d-DWT method was used. All of them correspond to a 4-5-3 band composite image.

3.1.1.

Landsat

The first plot of Fig. 5 is a representative result of the alterations in the spatial pattern caused by the lossy compression (BIFR in this case) mechanisms. The main alteration to the spatial pattern due to lossy compression is the reduction in data variability at high compression ratios (more than $100:1$), although at lower CRs (from $2.5:1$ to $50:1$) the variograms have very similar shapes and structural parameters (nugget, range and sill). Table 2 is a summary of the behavior of the total sill, which is the most representative parameter in this study, for two dates and two regions for the BIFR method. Comparing the sill between compressed and noncompressed images represents a measurement of the loss of detail and variability caused by a compression method for all Landsat spectral bands. Table 3 shows the same information in relation to the 3d-DWT compression method. As seen in Fig. 6, there are only small differences between the two methods for Landsat images with the highest CR. Therefore, the following tables and figures corresponding to Landsat images only show the results for the BIFR method. The next cases (CASI and MODIS) demonstrate that 3d-DWT is a more useful method for image series that are larger than Landsat image series.

Fig. 5

Variogram results of the BIFR compression method on a Landsat image from the Ripollès region on 02-07-2006, a multispectral view.

Table 2

Two examples (Ripollès region on 02-07-2006, and Penedès region on 11-9-2006) of the variation in the variogram parameter total sill at different compression ratios for the BIFR method.

	Band/CR	$1:1$	$2.5:1$	$5:1$	$10:1$	$20:1$	$100:1$	$200:1$	$400:1$
Ripollès02-07-2006	1-B	5.92	5.95	5.96	5.95	5.80	5.04	4.51	4.17
	2-G	5.30	5.34	5.35	5.33	5.23	4.58	4.13	3.81
	3-R	9.40	9.44	9.44	9.43	9.31	8.3	7.90	7.07
	4-nIR	92.96	93.30	93.34	93.19	90.93	82.33	74.19	69.57
	5-mIR1	26.31	26.40	26.43	26.38	25.77	22.63	21.28	18.85
	6-tIR	7.10	7.10	7.11	7.11	7.09	6.84	6.72	6.53
	7-mIR2	15.74	15.78	15.81	15.69	15.40	13.8	13.41	12.50
Penedès11-9-2006	1-B	18.01	18.51	18.57	18.53	18.06	15.68	15.68	15.68
	2-G	17.30	17.73	17.79	17.77	17.33	14.96	14.96	14.95
	3-R	25.94	26.37	26.43	26.31	25.68	22.21	22.21	22.21
	4-nIR	36.49	37.8	38.00	37.49	36.23	30.33	30.33	30.32
	5-mIR1	44.61	45.53	45.67	45.11	43.76	37.43	37.43	37.43
	6-tIR	3.53	3.54	3.55	3.53	3.53	3.35	3.35	3.36
	7-mIR2	31.54	32.06	32.17	31.9	31.15	27.24	27.23	27.24

Table 3

Two examples (Ripollès region on 02-07-2006, and Penedès region on 11-9-2006) of the variation in the variogram parameter total sill at different compression ratios for the 3d-DWT method.

	Band/CR	$1:1$	$2.5:1$	$5:1$	$10:1$	$20:1$	$100:1$	$200:1$	$400:1$
Ripollès02-07-2006	1-B	5.92	5.95	5.97	5.93	5.88	5.54	5.35	5.02
	2-G	5.30	5.33	5.34	5.33	5.30	5.05	4.89	4.61
	3-R	9.40	9.42	9.45	9.45	9.37	9.06	8.87	8.65
	4-nIR	92.96	92.80	92.83	92.19	90.44	82.75	76.49	70.17
	5-mIR1	26.31	26.37	26.40	26.04	25.33	22.47	20.42	18.34
	6-tIR	7.10	7.08	7.10	7.20	7.29	7.31	7.28	7.16
	7-mIR2	15.74	15.76	15.80	15.80	15.62	14.37	13.63	12.84
Penedès11-9-2006	1-B	18.01	18.45	18.52	18.59	18.38	16.81	16.37	15.69
	2-G	17.30	17.64	17.72	17.71	17.53	16.33	15.76	14.93
	3-R	25.94	26.26	26.35	26.06	25.34	22.47	20.73	19.15
	4-nIR	36.49	37.62	37.77	37.40	35.99	31.00	27.46	23.84
	5-mIR1	44.61	45.33	45.43	44.98	43.53	39.05	36.35	32.06
	6-tIR	3.53	3.52	3.56	3.82	4.28	5.18	5.73	6.44
	7-mIR2	31.54	31.93	32.03	32.34	32.33	30.53	29.44	27.97

Fig. 6

Landsat variogram comparison between BIFR and 3d-DWT compression in relation to the original pattern.

The sequenced plot in Fig. 5 shows that there are no significant differences in the variogram shape between bands, although there are differences in band range variances, and consequently in structural parameters depending on the region of the landscape.

Table 4 and Fig. 7 reproduce the same variogram spatial pattern in a temporal sequence for an example of the near infrared (4-nIR) band.

Table 4

Time comparison of variogram results (total sill) for BIFR compression of Landsat images from the Penedès region on the 1-B and nIR bands.

	Date	Band/CR	$1:1$	$2.5:1$	$5:1$	$10:1$	$20:1$	$100:1$	$200:1$	$400:1$
Penedès	13/04/2006	1-B	22.32	22.01	22.03	22.21	21.96	20.07	19.12	17.72
		4-nIR	38.46	38.48	38.48	38.23	36.88	31.36	28.10	24.63
	02/07/2006	1-B	20.51	20.12	20.12	20.29	20.21	18.68	18.37	17.42
		4-nIR	40.49	39.70	39.70	39.38	38.09	33.33	30.04	26.95
	19/08/2006	1-B	37.96	37.27	37.29	37.34	36.94	33.89	32.60	30.76
		4-nIR	54.20	53.35	53.38	53.18	51.74	45.71	42.01	37.39
	11/09/2006	1-B	18.01	18.45	18.52	18.59	18.38	16.81	16.37	15.69
		4-nIR	36.49	37.62	37.77	37.40	35.99	31.00	27.46	23.84

Fig. 7

Date comparison of variogram results for BIFR compression of Landsat images from the Penedès region on the near infrared (4-nIR) band.

Furthermore, lossy compression may produce different variance alterations in relation to the spatial direction. In order to analyze this possible anisotropy modification, variance pattern at different azimuth angles has been studied. Figure 8 shows that the lossy compression mechanisms used in this work maintain the expected directional patterns at all compression ratios; thus, this result proves that the previous variograms were always omnidirectional.

Fig. 8

Variogram behavior in relation to different azimuth directions. Landsat 02-07-2006 1-B band images.

3.1.2.

CASI

The main characteristic of CASI images in relation to Landsat images is their higher spectral resolution (much more data in the spectral domain). The higher spatial and radiometric resolution are also notable, but they are not the focus of the present study. It is therefore an interesting test bed for exploiting three-dimensional compression methods. As it is impossible to show the results for 72 bands for two regions and two dates, the following figures display representative results for selected bands for different regions and dates. For example, Table 5 shows the total sill parameter of the 3d-DWT method for four selected bands (10, 21, 34, and 58) in the $k4p$ scene (19-06-2007), while Fig. 9 plots variograms at different compression ratios for three selected bands. Figures 10 and 11 show the variogram analyses of the $c3p$ scene for the two regions applying the 3d-DWT compression method.

Table 5

Total sill parameters for the k4p scene in the Selva region and for the complete sequence of compression ratios at four representative spectral bands.

	Band/CR	$1:1$	$2.5:1$	$5:1$	$10:1$	$20:1$	$100:1$	$200:1$	$400:1$
SelvaK4p	10	14.07	14.07	14.07	13.82	13.91	12.25	10.37	8.71
	21	36.88	36.39	36.41	36.21	35.13	29.81	26.40	23.62
	34	40.92	41.18	41.20	40.97	40.26	37.71	35.80	31.88
	58	398.05	398.05	397.98	395.84	388.54	370.07	356.27	334.35

Fig. 9

Variogram results for the 3d-DWT compression method in the Selva region at selected CASI bands on 19-06-2007 ($k4p$ scene).

Fig. 10

Comparison of variogram results for the 3d-DWT compression method in the Selva region at selected CASI bands on 18-05-2007 ($c3p$ scene).

Fig. 11

Variogram results for the 3d-DWT compression method in the Vallès region at selected CASI bands.

Table 6 shows the total sill parameter for the other scene ($c3p$ on 18-05-2007) in the two study regions (Selva and Vallès) comparing the BIFR and 3d-DWT compression methods at selected CRs and bands (10, 21, 34). This table comparison demonstrates that spatial patterns have better fidelity with the 3d-DWT method than the BIFR method for CASI images. Figure 12 shows that 3d-DWT plots are closer to the noncompressed patterns, and that BIFR significantly alters the spatial patterns, especially for high compression ratios.

Table 6

Comparison between total sill parameters for the BIFR and 3d-DWT methods applied to the c3p scene in two regions at four representative compression ratios and three spectral bands.

$c3p$	Band/CR	$1:1$	$2.5:1$		$5:1$		$20:1$		$100:1$
			BIFR	3d	BIFR	3d	BIFR	3d	BIFR	3d
Selva	10	38.65	38.65	38.66	38.69	38.67	37.24	38.31	30.59	38.31
	21	84.98	84.99	84.98	78.85	85.01	78.85	83.59	65.12	83.59
	34	105.63	105.97	105.93	101.46	106.01	101.46	103.88	85.42	103.88
Vallès	10	65.88	65.88	65.87	66.05	65.84	63.90	65.61	55.97	62.49
	21	99.38	99.38	99.38	99.40	99.36	98.46	98.46	82.55	95.26
	34	139.98	142.59	139.99	143.54	139.99	139.11	138.73	124.53	134.03

Fig. 12

Comparison of variogram results of the 3d-DWT and BIFR compression methods in both regions at two selected CASI bands.

3.1.3.

MODIS

The MODIS test bed allows a much larger amount of data to be explored in the time domain than the Landsat images because the MODIS daily revisit period makes it possible to obtain a large cloud-free time series. These series are especially appropriate for analyzing the pattern alterations of three-dimensional compression methods versus two-dimensional ones.

The results for the MODIS images are similar to those obtained for CASI images, but both are different from the Landsat images. In these cases, the 3d-DWT method improves the quality of BIFR compressions. Table 7 shows that, at low compression ratios, 3d-DWT and BIFR maintain the spatial variability patterns of the original images, and that, at high compression ratios, BIFR loses quality compared to the 3d-DWT method. The plots in Figs. 13 and 14 (corresponding to the same region but to different spectral bands) show a reduction in the total sill parameter at increasing compression ratios for the 3d-DWT method, while Fig. 15 shows that 3d-DWT produces the pattern variability more faithfully than BIFR.

Table 7

Comparison between the BIFR and 3d-DWT methods of the total sill for four selected dates in the Penedès region for two spectral bands at representative compression ratios.

Penedès	CR	$1:1$	$2.5:1$		$20:1$		$70:1$
Date	Band		BIFR	3d	BIFR	3d	BIFR	3d
176	RED	20.04	20.06	20.02	19.10	19.10	12.50	19.45
	NIR	19.46	19.43	19.46	17.62	18.70	13.93	18.73
207	RED	20.23	20.27	20.21	20.20	20.20	11.67	19.70
	NIR	20.41	20.48	20.46	19.35	20.04	13.63	18.56
249	RED	19.77	19.86	19.92	19.90	19.90	11.05	18.52
	NIR	29.05	29.06	29.13	26.45	28.49	23.62	27.20
271	RED	16.43	16.42	16.43	15.83	15.83	10.13	13.97
	NIR	18.72	18.70	18.70	17.47	17.63	9.91	16.25

Fig. 13

Comparison of experimental variograms of the 3d-DWT compression method for the Penedès region on three selected dates applied to MODIS images in the red band.

Fig. 14

Variograms of the 3d-DWT compression method for the Penedès region on three selected dates for MODIS images in the near infrared band.

Fig. 15

Comparison of variograms for the BIFR and 3d-DWT methods for three selected dates in the Penedès region at representative compression ratios for MODIS images in the near infrared band.

In summary, these results quantify spatial pattern alterations for lossy compression images depending on compression ratios, compression methods, and series dimension through analyzing differences between structural variogram parameters (focusing on total sill) of compressed and noncompressed images. These results shows that, for compression ratios higher than $1:100$, the variogram is clearly degraded, reducing variability and thus, some capabilities or accuracy in related applications. This degradation could be partially solved, in some cases, by 3d-DWT, but this method needs a large amount of redundant data for significantly improve the BIFR method. These improvements are similar for time series and spectral series.

This significant alteration (CR higher than $1:100$) could be relevant in several geostatistical studies of remote sensing images, such as geostatistical applications referred in Sec. 1 of this work, and for example, in landscape scale classification based on variogram analysis.⁵² In this example, if someone uses compressed remote sensing images instead of original images, and this compression is not done with the appropriate parameters, it could significantly alter classification results and, consequently, its main goal of determining the optimal support size (i.e., optimal spatial resolution) for characterizing forest ecosystems.

3.2.

Performance Results

All the variogram analyses were processed using the IBM cluster of the research laboratory of the Computer Architecture and Operating Systems Department at the Universitat Autònoma de Barcelona. The IBM cluster is formed by 32 nodes, each with two Dual Core Intel(R) Xeon (R) 3.0 GHz processors with 12 Gbyte of RAM, and communicated with an integrated dual gigabit ethernet. Table 8 shows the average time of 14 independent executions using a different number of workers, while Fig. 16 shows the graphical representation and corresponding speedup evaluation.⁵³

Table 8

Execution time and speedup using n (0 to 24) workers.

Fig. 16

Execution times (a) and speedup (b) with a different number of workers. The red dotted line corresponds to the theoretical linear speedup, while the blue line and points correspond to the empirical data.

The proximity between the empirical speedup behavior and the theoretical linear speedup confirms that the parallel design and implementation provide a satisfactory solution for the computational problem at hand. These performance results demonstrate the validity of the proposed parallel solution, and the significant time reduction evidences the benefits of using distributed environments for processing large amounts of data, such as remotely sensed images. In fact, they provide an exhaustive test bed allowing running many executions with different parameters at several compression ratios in different regions on different dates, etc., and obtaining faster results and more agile comparisons.