k945
March 2, 2018, 9:27am
1
Hi,
I am ending up with
when trying to plot the following netcdf data
As far as I can tell I have chosen the appropriate projections.
using NetCDF, PyPlot, PyCall;
@pyimport seaborn as sns;
@pyimport cartopy.crs as ccrs;
@pyimport cartopy.feature as cfeature;
sns.set()
lats = ncread("T2.nc", "XLAT");
lons = ncread("T2.nc", "XLONG");
sst = ncread("T2.nc", "T2");
ax = axes(projection = ccrs.RotatedPole())
contourf(lons, lats, sst, transform = ccrs.PlateCarree())
ax[:set_xlim](-80, 10);
ax[:set_ylim](-50, 5);
Any ideas on how why this would be happening?
I think you need to convert the lats and lons to the projected coordinates of your map.
See here for example : https://github.com/Balinus/ClimateTools.jl/blob/master/src/mapping.jl
Around lines 60-62 for the projected transformation and lines 32-50 for the map definition. This is using Basemap though, I assume Cartopy is similar in their approach.
m = basemap[:Basemap](projection="cyl", llcrnrlat = -90, urcrnrlat = 90, llcrnrlon = 0, urcrnrlon = 360, resolution = "c")
lon2, lat2 = meshgrid(lon, lat)
x, y = m(lon2, lat2)
edit - Looked at your nc file and there is no time definition in your file.
I get the following error when trying to map your data. You might want to check you lons and lats matrix. I get a similar contour plot as you have shown.
cs = m[:contourf](x, y, sst)
WARNING: x coordinate not montonically increasing - contour plot
may not be what you expect. If it looks odd, your can either
adjust the map projection region to be consistent with your data, or
(if your data is on a global lat/lon grid) use the shiftgrid
function to adjust the data to be consistent with the map projection
edit - The lons and lats matrix are not consistent I think.
lons
223×152 Array{Float64,2}:
80.0567 80.3708 80.6809 80.9873 81.2899 81.5889 … 108.866 109.042 109.218 109.395 109.573 109.752
80.52 80.8341 81.1442 81.4505 81.753 82.0518 109.204 109.378 109.553 109.729 109.905 110.082
80.9855 81.2995 81.6096 81.9157 82.2181 82.5168 109.542 109.715 109.888 110.063 110.237 110.413
81.4532 81.7671 82.0771 82.3831 82.6853 82.9838 109.88 110.052 110.224 110.396 110.57 110.744
81.923 82.2369 82.5467 82.8526 83.1546 83.4529 110.218 110.389 110.559 110.731 110.903 111.075
82.3951 82.7089 83.0185 83.3242 83.6261 83.9241 … 110.557 110.726 110.895 111.065 111.236 111.407
82.8695 83.1831 83.4925 83.798 84.0996 84.3974 110.896 111.063 111.231 111.4 111.569 111.739
83.3461 83.6595 83.9688 84.2741 84.5754 84.8729 111.235 111.401 111.567 111.734 111.902 112.07
83.825 84.1382 84.4473 84.7523 85.0533 85.3505 111.574 111.739 111.904 112.069 112.236 112.402
84.3063 84.6193 84.9281 85.2327 85.5335 85.8303 111.914 112.077 112.24 112.405 112.569 112.735
⋮ ⋮ ⋱ ⋮
-148.565 -148.878 -149.187 -149.492 -149.793 -150.09 -176.314 -176.479 -176.644 -176.809 -176.976 -177.142
-148.086 -148.399 -148.709 -149.014 -149.315 -149.613 … -175.975 -176.141 -176.307 -176.474 -176.642 -176.81
-147.609 -147.923 -148.233 -148.538 -148.84 -149.137 -175.636 -175.803 -175.971 -176.14 -176.309 -176.478
-147.135 -147.449 -147.758 -148.064 -148.366 -148.664 -175.297 -175.466 -175.635 -175.805 -175.976 -176.147
-146.663 -146.977 -147.287 -147.593 -147.895 -148.193 -174.958 -175.129 -175.299 -175.471 -175.643 -175.815
-146.193 -146.507 -146.817 -147.123 -147.425 -147.724 -174.62 -174.792 -174.964 -175.136 -175.31 -175.484
-145.725 -146.039 -146.35 -146.656 -146.958 -147.257 … -174.282 -174.455 -174.628 -174.803 -174.977 -175.153
-145.26 -145.574 -145.884 -146.19 -146.493 -146.792 -173.944 -174.118 -174.293 -174.469 -174.645 -174.822
-144.797 -145.111 -145.421 -145.727 -146.03 -146.329 -173.606 -173.782 -173.958 -174.135 -174.313 -174.492
k945
March 2, 2018, 3:09pm
4
My understanding is that cylindrical equidistant is the equivalent to the rotatedpole projection in cartopy. But I will have another look in case I missed something in the cartopy documentation.
Thanks
Well, the projection is not important here. I was simply suggesting to transform your lat-lon matrix to projected coordinates.
With that being said, I think you problem arise from your lat-lon matrix values.
k945
March 2, 2018, 4:01pm
6
It looks like you were right, after looking through some ncl scripts, it turns out that the lons are adjusted; such that there are no negative values
lons[lons .< 0] = 360 + lons[lons .< 0]
Making that adjustment has fixed the issue
Thanks for the help
1 Like
Great!
What’s your colormap? I like it.
Is there a Julia native alternative to do this kind of mapping? Cartopy looks like a very useful tool but I would like to have the freedom of using Plots.jl
k945
March 11, 2018, 10:42am
10
Sorry I am not aware of an Julia alternative. With that said, plots.jl uses pyplots, so it with a little bit of work it should be possible to do this using plots.jl
GMT.jl can certainly do it, though not in one step because the nc file is not a simple grid but arrays of lon,lat,z
This example , although it still uses the the more hard-core syntax, may be helpful.