Matlab weather map with geoshow

10 views (last 30 days)
Jack Ts
Jack Ts on 16 Jan 2016
Commented: Jack Ts on 19 Jan 2016
Hello all! I have a certain problem creating a weather map. I have a netcdf file including all the data (Wind and coord variables etc) which I have assigned to variables and tried creating the map. I 'm including the code below
ncid = netcdf.open('lap_wrfout_d02_20130620.nc');
finfo = ncinfo('lap_wrfout_d02_20130620.nc');
disp(finfo);
dimNames = {finfo.Dimensions.Name};
dimMatch = strncmpi(dimNames,'x',1);
disp(finfo.Dimensions(dimMatch));
vartime = ncread('lap_wrfout_d02_20130620.nc', 'time');
vardataU10 = ncread('lap_wrfout_d02_20130620.nc', 'U10');
vardataV10 = ncread('lap_wrfout_d02_20130620.nc', 'V10');
meanU10 = mean(vardataU10, 3);
meanV10 = mean(vardataV10, 3);
ws = sqrt(meanU10.^2 + meanV10.^2);
XLONGdata = ncread('lap_wrfout_d02_20130620.nc', 'XLONG');
XLATdata = ncread('lap_wrfout_d02_20130620.nc', 'XLAT');
Now the problem is i cannot create the weather map. I 've tried that thus far:
[Xi, Yi] = meshgrid(linspace(10,45),linspace(30,50));
vq = griddata(double(XLONGdata),double(XLATdata),double(ws),Xi,Yi);
axesm('MapProjection','lambert','MapLatLimit',[min(min(double(XLATdata))) max(max(double(XLATdata)))],'MapLonLimit',[min(min(double(XLONGdata))) max(max(double(XLONGdata)))],'Frame','on','Grid','on', 'MeridianLabel', 'on', 'ParallelLabel', 'on');
imagesc(linspace(10,45),linspace(30,50),vq, 'AlphaData', 0.85);
colormap 'jet';
xlabel('Longitude');
ylabel('Latitude');
colorbar
h = geoshow(Yi,Xi,vq,'DisplayType','surface','FaceColor', 'white','FaceAlpha',0,'linewidth',2);
h.ZData = zeros(size(h.XData));
p = geoshow('landareas.shp','FaceColor', 'white','FaceAlpha',0,'linewidth',2);
The lambert map appears as a small black dot compared to the weather map data(which appears with correct colors etc.) Thanks in advance. Edit: would it be useful to use contour plot?

Accepted Answer

Chad Greene
Chad Greene on 18 Jan 2016
I see the problem--You're mixing regular cartesian coordinates with map coordinates.
When you call Mapping Toolbox plotting functions like plotm, geoshow, etc., they work by converting lats/lons to projected coordinates before plotting. As an example you can explore how linem works by typing
open linem
All the linem funciton does is convert lats,lons to x,y via mfwdtran, then plot with line. So
linem(lat,lon)
is equivalent to
[x,y] = mfwdtran(lat,lon);
line(x,y)
The range of x,y values produced by mfwdtran is usually somewhere between -1 to 1. However, in your code above you initialize map axes with axesm and then try to add gridded data with imagesc. You have to pick one or the other--either projected coordinates with axesm and geoshow or unprojected coordinates with imagesc.
It looks like your Xi and Yi values are actually lats and lons? You might want to change the names of Xi and Yi to loni and lati to make it more clear, then use pcolorm instead of imagesc:
[loni, lati] = meshgrid(linspace(10,45),linspace(30,50));
pcolorm(lati,loni,vq)
  1 Comment
Jack Ts
Jack Ts on 19 Jan 2016
Hello and thanks for your answer! I've actually been able to find a solution, using contourm instead of "turning" my coords into an "imagemap". The code is as follows :
figure;
colormap('jet');
axesm('MapProjection','lambert','MapLatLimit',[min(min(double(XLATdata))) max(max(double(XLATdata)))],'MapLonLimit',[min(min(double(XLONGdata))) max(max(double(XLONGdata)))],'Frame','on','Grid','on', 'MeridianLabel', 'on', 'ParallelLabel', 'on');
xlabel('Longitude');
ylabel('Latitude');
PlotMap=contourm(double(XLATdata), double(XLONGdata), double(ws),'Fill', 'on');
c=colorbar;
c.Label.String='Wind m/s';
p = geoshow('landareas.shp','FaceColor', 'white','FaceAlpha',0,'linewidth',2);
Once again thanks for your answer, it was very informative :)

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!