# "Not-so-good" NetCDF

## Data source

[ESA Permafrost Climate Change Initiative](https://data.ceda.ac.uk/neodc/esacci/permafrost/data/permafrost_extent/L4/area4/pp/v03.0), Permafrost extent for the Northern Hemisphere, v3.0

## QGIS

Adding data to the map as usual: the layer has incorrect georeference.

<figure><img src="/files/CfKRcL5YqP7xAurAhXGy" alt=""><figcaption><p>data view in QGIS</p></figcaption></figure>

## Python + XARRAY + GDAL

Let's try to study the dataset with Xarray

The full notebook download is below. To run it, [set up conda environment](/working-with-scientific-data-formats-in-gis/intro-notes.md#configuring-environment-for-jupyter-notebook), savethe notebook to your working directory, and run Jupyter Notebook from this directory.&#x20;

{% file src="/files/8McIHnfW09URjwrIv5Fu" %}

### Read and view data

{% code overflow="wrap" %}

```python
import xarray as xr

ds = xr.open_dataset('ESACCI-PERMAFROST-L4-PFR-MODISLST_CRYOGRID-AREA4_PP-2019-fv03.0.nc')
```

{% endcode %}

<figure><img src="/files/XTopPfTy27USAKZK4qWM" alt=""><figcaption><p>Dataset has 2 variables, one of which is crs, but x and y coordinates are presented as 1-d arrays</p></figcaption></figure>

<pre class="language-python"><code class="lang-python"><strong>prf = ds.PFR[0,:,:]
</strong>plt.imshow(prf, cmap='Spectral', vmin=0, vmax=100)
plt.colorbar()
</code></pre>

<figure><img src="/files/019jc2rkfoBs2HoKvaJg" alt=""><figcaption><p>data view: the rotation is correct</p></figcaption></figure>

To export data as GeoTIFF some info georeferencing is needed. For GDAL it is a projection (crs reference) and [geotransform ](https://gdal.org/tutorials/geotransforms_tut.html)parameters (upper left x and y coordinates and pixel size). Some of this info is stored in the `polar_stereographic` variable.

<figure><img src="/files/ogQgCJ8R8cM9FihxWf1F" alt=""><figcaption></figcaption></figure>

`crs_wkt` is in GDAL-compatible format, we can use it. But `GeoTransform` is not consistent with dataset upper left pixel coordinates:

{% code overflow="wrap" %}

```python
print(ds.polar_stereographic.GeoTransform)
print(ds.x.values.min(), ds.x.values[1] - ds.x.values[0], 0, ds.y.values.max(), 0, ds.y.values[1] - ds.y.values[0])
> -8679599.425969256 926.6254331383326 0 4120958.560533902 0 -926.6254331383326 
> -6111475.222394747 926.6254331385717 0 4114895.094696621 0 -926.625433138106
```

{% endcode %}

### Save to GeoTIFF

{% code overflow="wrap" %}

```python
def save_file_singleband(input_array, out_file, projection, geotransform):
    
    """Save numpy array to singleband GTiff file"""
    
    driver = gdal.GetDriverByName('GTiff')
    
    dataset = driver.Create(out_file,
                        input_array.shape[1],
                        input_array.shape[0],
                        eType=gdal.GDT_Float32)
    
    #set projection
    dataset.SetProjection(projection)
    # set geotransform
    dataset.SetGeoTransform(geotransform)

    outband = dataset.GetRasterBand(1)
    outband.WriteArray(input_array)
    
    dataset = None
```

{% endcode %}

```python
save_file_singleband(
    prf,
    'permafrost_proj_wkt.tif',
    ds.polar_stereographic.crs_wkt,
    [-6111475.22239475,926.62543314,0,4114895.09469662,0,-926.62543314]
)
```

<figure><img src="/files/ZKJ0JYaTETDLvLiQEdBM" alt=""><figcaption><p>Output file in QGIS map view</p></figcaption></figure>


---

# Agent Instructions: 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:

```
GET https://guides.geospatial.bas.ac.uk/working-with-scientific-data-formats-in-gis/not-so-good-netcdf.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
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.