Analyzing 32 Billion Points of Chinese PM2.5 Readings

As the world continues to advance technologically, our natural environment silently bears the cost. Among the most critical indicators of this deterioration is the air we breathe. Many countries today face dangerously high pollution levels, posing serious health risks to millions.

This growing problem calls for deeper understanding — because identifying a problem is the first step toward solving it. In this project, we explore a large-scale dataset of Chinese PM2.5 readings, which measure the concentration of fine particulate matter (particles smaller than 2.5 micrometers) suspended in the air. The higher the PM2.5 concentration, the more harmful the air is to breathe.

We obtained the dataset from **Zenodo, spanning the years 2021 to 2024**. Our objective was to analyze and visualize PM2.5 patterns across China — to understand spatial and temporal pollution trends and lay the groundwork for forecasting models that could support better environmental planning and policy.

This post focuses on the Exploratory Data Analysis (EDA) phase of that process.

Understanding the Dataset

According to the **Zenodo documentation, the dataset provides daily PM2.5 readings across China at a spatial resolution of approximately 1 km × 1 km — around 22 million data points per day, and more than 32 billion** measurements in total for the four years combined.

After downloading and extracting the ZIP archives (each containing multiple .nc file per day), our first step was to verify dataset completeness — ensuring that we had all daily files across the intended years. Below is the picture that ensures the completeness of the files w.r.t. years.

Verification for completeness

Once completeness was confirmed, the next step was to inspect the structure and type of data within each file. Specifically, we needed to confirm that all files followed a uniform format and contained consistent variables. For this, we checked that every .nc file:

• Contained the same total number of data points (i.e., identical grid size and spatial resolution).
• Included the same set of variables and features.

Below are the main structural characteristics of the dataset:

1. Each file name encodes the date in the format year_month_day.
2. The core features are Latitude, Longitude, and PM2.5.
3. Latitude and longitude values are spaced at 0.01°, representing roughly 1 km² grid cells.
4. PM2.5 readings are measured in µg/m³ (micrograms per cubic meter).

After confirming these structural details, we proceeded to analyze the distribution of missing values (NaNs) across the dataset. This step was critical to determine whether the missing data followed a consistent spatial pattern across files. If NaNs appeared at the same grid locations in all files, it would likely indicate areas outside China’s boundaries or regions not covered by sensors. However, if NaNs varied between files, it would suggest incomplete measurements that might require imputation or exclusion in later steps.

Because the dataset consists of multi-dimensional NetCDF files, we utilized Xarray, a library specifically designed for labeled, multi-dimensional data. Xarray operates on top of Dask, which enables lazy evaluation and parallel processing. Unlike traditional libraries such as Pandas that load data entirely into memory, Xarray with Dask processes data in manageable chunks — allowing us to efficiently handle massive datasets like this one without exhausting system resources.

All .nc files were read lazily using Xarray, allowing efficient processing of large datasets without exhausting memory resources. For each file, key information such as dimensions, variables, and missing value counts was extracted and stored in a structured DataFrame. This enabled us to verify that all files shared the same structure, including identical features, dimensions, and the number of daily entries. We then calculated the average number of missing values per file and per variable and summarized these statistics in CSV files for further validation. To determine whether missing values occurred at the same spatial indices, we generated missing-value masks for each file. However, since each mask contained around 22 million grid points, direct comparison was computationally expensive. To address this efficiently, we utilized the Hashlib library to create compact hash representations of these masks. By comparing the hashes, we confirmed that all files shared an identical missing-value pattern. Therefore, since both the structure and missing-value distributions were consistent across all files, no imputation was necessary, and the dataset was ready for further exploratory data analysis (EDA). Below is the CSV output we got after analyzing all of the files.

DATA COMPLETNESS RESULTS

Now lets jump to EDA!

How do PM2.5 concentrations vary spatially across different provinces in China and temporally across seasons or months?

To analyze how PM2.5 concentrations vary spatially across Chinese provinces and temporally across seasons or months, we implemented an efficient geospatial extraction pipeline using GeoPandas, Xarray, and NumPy with memory-safe batch processing.

1. Province Sampling: We began by loading China’s provincial boundaries from the GADM GeoJSON file. For each province, we generated 50 random points located strictly within its boundaries using geometric sampling. These serve as representative spatial locations for PM2.5 measurement.
2. Grid Validation: Since PM2.5 data were stored in gridded .nc files, we mapped each sampled point to its nearest latitude–longitude grid cell using vectorized nearest-neighbor matching. Points beyond a 20 km distance threshold from any valid grid cell were excluded, ensuring spatial validity.
3. Province–Grid Index Mapping: The validated grid indices for each province were stored in a dictionary (and serialized as valid_location_indices.pkl) for efficient repeated access without reprocessing.
4. Daily PM2.5 Extraction (Memory-Safe): Using an optimized function with direct NumPy indexing, we extracted the daily average PM2.5 values for all valid province locations from each .nc file spanning 2021–2024. To handle the massive dataset efficiently, the files were processed in batches of 50 with parallel workers, while intermediate results were periodically stored in Parquet format. Additionally, garbage collection and memory monitoring were employed throughout the process to prevent memory overflow and ensure stable execution.
5. Final Aggregation and Storage: The aggregated daily averages were compiled into pm25_daily_provinces.parquet and .csv, giving a continuous time series of daily PM2.5 per province.
6. The processed data was then visualized using Plotly, offering an intuitive view of the spatial and temporal variations in PM2.5 concentrations across provinces. These visualizations help reveal how the air quality in a specific province changes over time.

Spatiotemporal analysis

Where are the major pollution hotspots in China?

All NetCDF files were processed lazily with Dask for memory efficiency, where each file was read, downsampled, and flattened only over valid land regions defined by a NaN mask.

To reflect the natural diffusion of smog — similar to any gaseous dispersion — we adopted a 500km × 500 km spatial resolution instead of the raw 1 km grid. This subsampling not only reduced computational overhead but also provided a more meaningful spatial scale for atmospheric analysis. This all means we are assuming the conditions of smog are almost same in 500km × 500 km area.

Once all valid data points were stacked, we computed the mean PM₂.₅ concentration across time, identified the Top 10 pollution hotspots, and generated summary statistics. The resulting hotspots highlight the regions with persistently high particulate matter concentrations, revealing key spatial patterns of air pollution across the country. The top 10 places with highest PM₂.₅ were then plotted.

From the visualizations, we can observe that the top 10 regions are scattered across China, showing no single dominant cluster. An interesting pattern, however, is that areas located closer to Pakistan or India tend to exhibit poorer air quality. This aligns with broader regional trends — Pakistan consistently ranks among the countries with the highest PM2.5 concentrations, with cities like Lahore frequently appearing among the most polluted in the world.

Left: Top 10 locations with highest PM2.5 readings. Right: Zoomed in version

Which regions have improved and which regions have worsened their air quality over time?

This question helps us pinpoint which regions have improved or deteriorated in air quality over time. Identifying these locations allows us to focus on the underlying factors that might have contributed to their respective trends — whether policy changes, industrial shifts, or climatic influences.

To ensure fair comparison across years, we applied the same data filtering and preprocessing strategy: removing missing (NaN) values and using smart spatial down sampling to handle vast datasets efficiently.

We then quantified the change for each region using a simple formulation:

Change = PM2.52021−PM2.52024​

• Positive change→ indicate improvement (air got cleaner).
• Negative change→ indicate worsening (air got more polluted).

Finally, we ranked all valid grid locations by this change metric and highlighted the Top 10 most improved and Top 10 most worsened regions for visualization. It can be seen that the location nearer to Pakistan and India are improving (which is a good news to what we saw in our earlier analysis).

Top 10 locations that have improved or degraded within the period 2021–2024

What is the spatial and temporal distribution of extreme pollution events over 4 years? How long do they last and are they becoming more frequent?

We applied the same smart sampling strategy again, and then identified all extreme pollution events where PM₂.₅ levels exceeded 150 µg/m³ over the four-year period (2021–2024). Each year’s data was processed in parallel Dask batches, and only the relevant extreme readings were extracted and stored in optimized Parquet files for efficiency.

This approach significantly reduced computation time and memory usage, enabling us to pinpoint when and where these extreme events occurred, how frequently they appeared, and how long they persisted. By aggregating these events across space and time, we could analyze whether such high-pollution episodes are becoming more frequent or prolonged in particular regions.

From the results, it can be observed that locations closer to Pakistan and India tend to exhibit the most severe pollution levels, which aligns with earlier findings. Overall, however, the extreme-event locations appear to be fairly evenly distributed across China.

Left: Shows the location on map. Right: shows the number of events in each year

How do PM2.5 concentrations differ between major urban centers and rural areas in China over the time. Is the urban-rural pollution gap widening or narrowing?

We started by defining the latitude and longitude of 10 urban and 10 rural locations across China. Using these coordinates, we lazily searched through all NetCDF PM2.5 files from 2021 to 2024, extracting values efficiently without loading entire datasets into memory. The extracted data was then compiled into a DataFrame, cleaned, and used to calculate yearly mean PM2.5 levels for urban and rural areas. Finally, we compared these means to quantify the urban-rural pollution gap and analyzed its trend over time, highlighting whether the gap is widening, narrowing, or remaining stable.

From the graphs, it is clear that urban areas generally have higher PM2.5 levels than rural areas, likely due to factors such as fewer trees, higher numbers of automobiles, and industrial activity. Overall, the difference between urban and rural areas has remained relatively consistent on average.

Rural and urban comparisons

Summary:

We performed exploratory data analysis (EDA) on a large PM2.5 dataset using Xarray with lazy evaluation, enabling efficient handling of massive NetCDF files without loading everything into memory. Along the way, we made smart, computation-saving decisions that optimized both time and resources. Our analysis also focused on interpreting the data in the context of specific questions, with insights clearly communicated and supported through visualizations.