Search code examples
pythonnetcdfmatplotlib-basemapcartopy

Plotting .netCDF data correctly, impossible in Python (Basemap, Cartopy)

I have been using .netCDF files (satellite data) for a few weeks now, and in general, I never had issues plotting them. I've been experimenting with sea ice concentration data but I can't get them right. As far as I can tell, the problem has to do something with the longitude/latitude being in meters. Nothing works right. Data seem to be rotated or falsely projected. The EPSG is 3411. I usually plot with cartopy but I decided to use Basemap also. None of them got the job done with my code. File can be found here.

import os
import xarray as xr
import cartopy.crs as ccrs 
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
from pyproj import Transformer

data = xr.open_dataset('df.nc')
sea_ice = data.F17_ICECON; sea_ice

lons = data.x
lats = data.y
lon, lat = np.meshgrid(lons, lats)

# Also experimented with this but without any results
transformer = Transformer.from_crs("EPSG:3411", "EPSG:4326")
lon_T, lat_T = transformer.transform(lon, lat)

# With Basemap
m = Basemap(width=6000000,height=6000000,
            resolution='l',projection='stere',lat_0=90,lon_0=0, epsg=3411)

m.pcolormesh(lon, lat, sea_ice[2])
m.drawcoastlines()
m.fillcontinents(color='grey')

# With Cartopy
day = 4
ax = plt.axes(projection=ccrs.PlateCarree())

plt.contourf(lons, lats, sea_ice[day],
             transform=ccrs.PlateCarree())

ax.coastlines()

plt.show()

This is the image produced with Basemap and this is with Cartopy. Actual data should look like this.

Solution

  • Just use Cartopy with appropriate projection and get the plot. The projection I use is for plotting/visualization purposes only, as it is not exact equal EPSG:3411.

    import xarray as xr
import cartopy.crs as ccrs 
import matplotlib.pyplot as plt
import numpy as np

data = xr.open_dataset('df.nc')
sea_ice = data.F17_ICECON

lons = data.x
lats = data.y
lon, lat = np.meshgrid(lons, lats)

fig = plt.figure(111, figsize=[8,8])

# Use this projection instead of the one that fails
useproj = ccrs.NorthPolarStereo(central_longitude=-45,true_scale_latitude=70)
ax = plt.subplot(projection=useproj)
day = 4
ax.contourf(lons, lats, sea_ice[day])
ax.coastlines(alpha=0.5, lw=0.3)
gl = ax.gridlines(draw_labels=True)

plt.show()

    npolarplot