plotting netcdf data on a map in matplotlib

jmssnyder · March 30, 2017, 8:16pm

I am trying to plot data from a netcdf file onto a map through matplotlib.
I initially tried to use basemap but then someone suggested that I use
cartopy based on the website:

http://scitools.org.uk/cartopy.

Either way I am able to generate the map fine but am having difficulty
plotting the data successfully. For instance using basemap I wind up
getting a weird plot of data that does not make sense. Also how to I go
about plotting the netcdf data around the dateline? Attached is the ncdump
output (gridoutput.txt) for the netcdf file in terms of the grid and also
the script that I am using to try to plot the data (netcdfplot.py).
Can anyone give me some suggestions on how to go about fixing the errors
that I have in the script?

-Jason
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/6e426fcf/attachment-0001.html>
-------------- next part --------------
netcdf \2017032800 {
dimensions:
  latitude = 159 ;
  longitude = 361 ;
  time = UNLIMITED ; // (1 currently)
variables:
  double latitude(latitude) ;
    latitude:units = "degrees_north" ;
    latitude:long_name = "latitude" ;
  double longitude(longitude) ;
    longitude:units = "degrees_east" ;
    longitude:long_name = "longitude" ;
  double time(time) ;
    time:units = "seconds since 1970-01-01 00:00:00.0 0:00" ;
    time:long_name = "verification time generated by wgrib2 function verftime()" ;
    time:reference_time = 1490659200. ;
    time:reference_time_type = 1 ;
    time:reference_date = "2017.03.28 00:00:00 UTC" ;
    time:reference_time_description = "analyses, reference date is fixed" ;
    time:time_step_setting = "auto" ;
    time:time_step = 0. ;
  float HTSGW_surface(time, latitude, longitude) ;
    HTSGW_surface:_FillValue = 9.999e+20f ;
    HTSGW_surface:short_name = "HTSGW_surface" ;
    HTSGW_surface:long_name = "Significant Height of Combined Wind Waves and Swell" ;
    HTSGW_surface:level = "surface" ;
    HTSGW_surface:units = "m" ;

// global attributes:
    :Conventions = "COARDS" ;
    :History = "created by wgrib2" ;
    :GRIB2_grid_template = 0 ;
data:

latitude = -78, -77, -76, -75, -74, -73, -72, -71, -70, -69, -68, -67, -66,
    -65, -64, -63, -62, -61, -60, -59, -58, -57, -56, -55, -54, -53, -52,
    -51, -50, -49, -48, -47, -46, -45, -44, -43, -42, -41, -40, -39, -38,
    -37, -36, -35, -34, -33, -32, -31, -30, -29, -28, -27, -26, -25, -24,
    -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9,
    -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
    13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
    31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
    49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66,
    67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80 ;

longitude = 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192,
    193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206,
    207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,
    221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234,
    235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248,
    249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262,
    263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276,
    277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290,
    291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304,
    305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318,
    319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332,
    333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346,
    347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 0, 1, 2,
    3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
    23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
    41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58,
    59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
    77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94,
    95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109,
    110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123,
    124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137,
    138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151,
    152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165,
    166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180 ;

time = 1490659200 ;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: erroneous_chart.png
Type: image/png
Size: 17376 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/6e426fcf/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: netcdfplot.py
Type: text/x-python
Size: 909 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/6e426fcf/attachment-0001.py>

jklymak1 · March 30, 2017, 8:38pm

I?m not an expert on cartopy.

1) can you plot this data w/o using basemap?
2) why do you have the ?squeeze? command in there?
3) I?m not clear what the problem is with your image from the snippet you sent. Is it just an issue with the longitude wrapping? Some projections demand longitude must be between -180 and 180. i.e. lon[lon>=180]=lon[lon>=180]-360. You may also need to re-order your longitude and data so that longitude is monotonically increasing.
4) I recommend xarray when dealing w/ netcdf files, though maybe that is overkill for your problem.

A minimal working example with a link to the bad data always helps.

I?m sure with a bit of experimentation and use of google you can get this working.

i.e.
http://stackoverflow.com/questions/13856123/setting-up-a-map-which-crosses-the-dateline-in-cartopy

Cheers, Jody

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/ac932038/attachment.html>

jmssnyder · March 30, 2017, 8:49pm

The main issue I have is that the data does not plot properly. There
evidently is an issue of matching coordinates in the data versus the map.
I am sending the netcdf file along with the script. I am not sure how to
plot a map with data without using basemap since this is what provides the
continental outlines. If you know of any way I can plot this data using
matplotlib successful let me know.

-Jason

···

On Thu, Mar 30, 2017 at 8:38 PM, Jody Klymak <jklymak at uvic.ca> wrote:

I?m not an expert on cartopy.

1) can you plot this data w/o using basemap?
2) why do you have the ?squeeze? command in there?
3) I?m not clear what the problem is with your image from the snippet you
sent. Is it just an issue with the longitude wrapping? Some projections
demand longitude must be between -180 and 180. i.e.
lon[lon>=180]=lon[lon>=180]-360. You may also need to re-order your
longitude and data so that longitude is monotonically increasing.
4) I recommend xarray when dealing w/ netcdf files, though maybe that is
overkill for your problem.

A minimal working example with a link to the bad data always helps.

I?m sure with a bit of experimentation and use of google you can get this
working.

i.e.
http://stackoverflow.com/questions/13856123/setting-up-
a-map-which-crosses-the-dateline-in-cartopy

Cheers, Jody

--
Jason Snyder PhD
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/300da166/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2017032800.nc.tar.gz
Type: application/x-gzip
Size: 30496 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/300da166/attachment-0001.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testnetcdf.py
Type: text/x-python
Size: 28 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/300da166/attachment-0001.py>

rmay31 · March 30, 2017, 9:36pm

Here's a version that works for me with cartopy. It relies on np.roll(),
which just got added in numpy 1.12, to correct the fact that the longitudes
are arranged 180-360, 0-180, which seems to be causing problems. The other
easy way to solve that is to plot the data in chunks.

import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import numpy as np

ncfile='/Users/rmay/Downloads/2017032800.nc'
fh = Dataset(ncfile)

lats = fh.variables['latitude'][:]
lons = fh.variables['longitude'][:]
HTSGW_var = fh.variables['HTSGW_surface']

roll_to = -lons.argmin()
lons = np.roll(lons, roll_to)
data = np.roll(HTSGW_var[:].squeeze(), roll_to, axis=-1)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_global()
lon, lat = np.meshgrid(lons, lats)
cs = ax.pcolor(lon, lat, data)
# If roll not available
#cs = ax.pcolor(lon[:, :180], lat[:, :180], data[:, :180])
#cs = ax.pcolor(lon[:, 180:], lat[:, 180:], data[:, 180:])

# draw coastlines.
ax.coastlines()
ax.add_feature(cfeat.LAND)
plt.show()

···

On Thu, Mar 30, 2017 at 2:49 PM, Jason Snyder <jmssnyder at ucdavis.edu> wrote:

The main issue I have is that the data does not plot properly. There
evidently is an issue of matching coordinates in the data versus the map.
I am sending the netcdf file along with the script. I am not sure how to
plot a map with data without using basemap since this is what provides the
continental outlines. If you know of any way I can plot this data using
matplotlib successful let me know.

-Jason

On Thu, Mar 30, 2017 at 8:38 PM, Jody Klymak <jklymak at uvic.ca> wrote:

I?m not an expert on cartopy.

1) can you plot this data w/o using basemap?
2) why do you have the ?squeeze? command in there?
3) I?m not clear what the problem is with your image from the snippet you
sent. Is it just an issue with the longitude wrapping? Some projections
demand longitude must be between -180 and 180. i.e.
lon[lon>=180]=lon[lon>=180]-360. You may also need to re-order your
longitude and data so that longitude is monotonically increasing.
4) I recommend xarray when dealing w/ netcdf files, though maybe that is
overkill for your problem.

A minimal working example with a link to the bad data always helps.

I?m sure with a bit of experimentation and use of google you can get this
working.

i.e.
http://stackoverflow.com/questions/13856123/setting-up-a-
map-which-crosses-the-dateline-in-cartopy

Cheers, Jody

--
Jason Snyder PhD

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users at python.org
Matplotlib-users Info Page

--
Ryan May
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170330/2556d05e/attachment.html>