> For the complete documentation index, see [llms.txt](https://jonas.gitbook.io/seismo_rain/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://jonas.gitbook.io/seismo_rain/rainfall-data/extracting-rainfall-data.md).

# Extracting rainfall data

Due to the remote location of the seismic network and the lack of nearby weather stations with rain gauges, as well as the highly dynamic nature of rainstorms on a local scale (less than a kilometre), we will use a satellite-derived precipitation reanalysis model to extract the rainfall data of interest.

## [**CPC Morphing Technique (CMORPH) High Resolution Global Precipitation Estimates, Version 1** ](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C00948)

The bias-corrected CMORPH dataset provides high-resolution global satellite precipitation estimates on an 8 km x 8 km grid, updated at 30-minute intervals. These estimates are derived from integrating passive microwave measurements and corrected using rain gauges and other numerical models.

To construct a robust result, this is the highest temporal and spatial resolution of rainfall data that we can achieve in this instance.

## Opening precipitation NetCDF files

To obtain precipitation estimates from the CMORPH\_V1.0 dataset for the period from April 19, 2016, UTC 22:00 to April 20, 2016, UTC 02:00, four netCDF (.nc) files are needed.&#x20;

These files can be downloaded from the following link: [CMORPH\_V1.0\_ADJ\_8km-30min Dataset](https://www.ncei.noaa.gov/data/cmorph-high-resolution-global-precipitation-estimates/access/30min/8km/). Make sure to change the file path to your local directory when opening it in `xarray`.

```python
import xarray as xr

#Precipitation data 2016-04-19 UTC 22:00 - 23:00
prep22 = xr.open_dataset("/Users/laichungheijonas/Downloads/CMORPH_V1.0_ADJ_8km-30min_2016041922.nc")

#Precipitation data 2016-04-19 UTC 23:00 - 00:00(+1)
prep23 = xr.open_dataset("/Users/laichungheijonas/Downloads/CMORPH_V1.0_ADJ_8km-30min_2016041923.nc")

#Precipitation data 2016-04-20 UTC 00:00 - 01:00
prep00 = xr.open_dataset("/Users/laichungheijonas/Downloads/CMORPH_V1.0_ADJ_8km-30min_2016042000.nc")

#Precipitation data 2016-04-20 UTC 01:00 - 02:00
prep01 = xr.open_dataset("/Users/laichungheijonas/Downloads/CMORPH_V1.0_ADJ_8km-30min_2016042001.nc")
```

## Concatenate all netCDF files&#x20;

Using `xr.concat`, we can concatenate all the netCDF files along with the time dimension.

```python
prep = xr.concat([prep22,prep23,prep00,prep01],dim='time')
print(prep)
```

Upon executing `print(prep)`, you should observe the following output. Given the 4-hour interval, there should be a total of 8 timesteps.

```
<xarray.Dataset> Size: 262MB
Dimensions:      (time: 8, nv: 2, lat: 1649, lon: 4948)
Coordinates:
  * time         (time) datetime64[ns] 64B 2016-04-19T22:00:00 ... 2016-04-20...
  * lat          (lat) float64 13kB -59.96 -59.89 -59.82 ... 59.82 59.89 59.96
  * lon          (lon) float64 40kB 0.03638 0.1091 0.1819 ... 359.8 359.9 360.0
Dimensions without coordinates: nv
Data variables:
    time_bounds  (time, nv) datetime64[ns] 128B 2016-04-19T22:00:00 ... 2016-...
    lat_bounds   (time, lat, nv) float64 211kB -60.0 -59.93 ... 59.93 60.0
    lon_bounds   (time, lon, nv) float64 633kB 0.0 0.07276 ... 359.9 360.0
    cmorph       (time, lat, lon) float32 261MB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes: (12/57)
    ncei_template_version:      NCEI_NetCDF_Grid_template_V2.0
    title:                      NOAA Climate Data Record (CDR) of CPC Morphin...
    keywords:                   Precipitation, Satellite, High-Resolution, Gl...
    summary:                    The CMORPH CDR is a reprocessed and bias-corr...
    references:                 Xie, P., et al. (2017), Reprocessed, Bias-Cor...
    Conventions:                CF-1.6, ACDD-1.3
    ...                         ...
    geospatial_lat_resolution:  0.072771376
    geospatial_lat_units:       degrees_north
    geospatial_lon_min:         0.0
    geospatial_lon_max:         360.0
    geospatial_lon_resolution:  0.072756669
    geospatial_lon_units:       degrees_east
```

## Extracting the rainfall data for each seismic station

To obtain the coordinates of each seismic station, we can utilize the inventory `inv` previously defined. The inventory can be divided into different levels as follows:

* `inv[0]`: Seismic Network Details
* `inv[0][0]`: Seismic Station Details
* `inv[0][0][0]`: Channel Details

```python
for x in range(4):
    print(f'{inv[0][x].code}: latitude: {inv[0][x].latitude}, longitude: {inv[0][x].longitude}')
```

Output:

```
1517: latitude: 36.796766, longitude: -98.105601
199: latitude: 36.796824, longitude: -98.020063
465: latitude: 36.7962, longitude: -97.92993
763: latitude: 36.796074, longitude: -97.839733
```

Subsequently, we can use the `.sel` command to extract the rainfall data for the specific grid cell. It's important to note that we specify `method='nearest'` to identify the nearest rainfall data grid cell to the seismic station.

`prep['cmorph']` is the variable for rainfall intensity in mm/hour/

```python
prep_1517 = prep['cmorph'].sel(lat=36.796766, lon=360-98.105601 , method='nearest')
prep_199 = prep['cmorph'].sel(lat=36.796824, lon=360-98.020063, method='nearest') 
prep_465 = prep['cmorph'].sel(lat=36.7962, lon=360-97.92993, method='nearest')
prep_763 = prep['cmorph'].sel(lat=36.796074, lon=360-97.839733, method='nearest')
```

## Plotting the rainfall intensity at all seismic stations

Since rainfall intensity is measured per half-hourly interval, the data point at 22:00 represents the aggregated rainfall intensity during the period from 22:00 to 22:30, and so forth for subsequent intervals.

```python
plt.figure(figsize=(15,3))

plt.plot(prep_1517.time,prep_1517,linewidth=3,label='2A.1517')
plt.plot(prep_199.time,prep_199,linewidth=3,label='2A.199')
plt.plot(prep_465.time,prep_465,linewidth=3,label='2A.465')
plt.plot(prep_763.time,prep_763,linewidth=3,label='2A.763')

plt.legend()
plt.ylabel("Precipitation (mm/hour)")
plt.xlabel("Time")
plt.show()
```

<figure><img src="/files/AeksIcN2hpC196UyOkCG" alt=""><figcaption><p>Rainfall Intensity at all four seismic station during 19 April 2016 UTC 22:00 - 20 April 2016 UTC 02:00</p></figcaption></figure>


---

# Agent Instructions
This documentation is published with GitBook. GitBook is the documentation platform designed so that both humans and AI agents can read, navigate, and reason over technical content effectively. Learn more at gitbook.com.

## Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter, and the optional `goal` query parameter:

```
GET https://jonas.gitbook.io/seismo_rain/rainfall-data/extracting-rainfall-data.md?ask=<question>&goal=<endgoal>
```

`ask` is the immediate question: it should be specific, self-contained, and written in natural language.
`goal` is optional and describes the broader end goal you are ultimately trying to accomplish on behalf of the user. GitBook uses it to tailor the answer towards what is most useful for that goal.

The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.
