Interpolate scattered latitude and longitude (climate) data and plot

40 views (last 30 days)
Hello everyone,
I have a table with 3 columns: latitude, longitude, and climate information. I want to plot these points on my map and have a color bar.
I want to interpolate them to achieve a zoning map.
I have a figure of the country as you can see below:
Using this code:
borders('Iran Islamic Republic of') % borders is a 3rd party function from matlab file exchange by Chad A. Greene
grid on
ax.GridLineStyle = '-';
% Modify the X and Y ticks positions
xticks(44:.5:65);
yticks(25:.5:40);
xlabel('Longitude')
ylabel('Latitude')
xtickangle(90)
set(gca,'FontSize',5)
set(hTitle,'FontSize',16)
xlabel('Longitude','FontSize',10,'FontWeight','bold','Color','r')
ylabel('Latitude','FontSize',10,'FontWeight','bold','Color','r')
I was tried to plot the value of latitude and longitude on a map using this code:
hold on
geoshow(SPI_for_certain_date.lat, SPI_for_certain_date.lon, SPI_for_certain_date.SPI_1month, 'displaytype', 'texturemap');
colorbar('vert','FontSize',12);
But it does not work and shows me this warning:
Warning: Error creating or updating Surface
Error in value of property ZData
Array is wrong shape or size
Warning: Error creating or updating Surface
Error in value of property ZData
Array is wrong shape or size
And so here my figure after that:
Now I don't know what should I do. The color bar understands the range of my data accurately but unfortunately data not plotted on my map. To check if everything is okay I plot just latitude and longitude on my map (without corresponding rainfall data) using pcolor and it works and points presented on a map.
Actually I want to do something like this code that I found it from Google (Internet):
% some sample data, located in europe
data = [randi([-10 30],10,1) randi([30 50],10,1) rand(10,1)];
lat = data(:,1);
lon = data(:,2);
temp = data(:,3); % : not ;
lon0 = min(lon) ; lon1 = max(lon) ;
lat0 = min(lat) ; lat1 = max(lat) ;
N = 100 ;
x = linspace(lon0,lon1,N) ;
y = linspace(lat0,lat1,N) ; % lat0, not lon1
[X,Y] = meshgrid(x,y) ;
F = scatteredInterpolant(lon,lat,temp) ;
Z = F(X,Y) ;
worldmap('World')
load coastlines
plotm(coastlat,coastlon)
geoshow(X,Y,Z,'DisplayType', 'texturemap') % options give it a nice colormap for temperature - but you might want to change them to something more of your taste
Here is world map I want to plot my file that I attahced on my map that I mentiond above.
As the original file (latitude longitude and values) is very large for attaching here, I cut part of it and so I attache an example file with a small size here.
Thank you so much in advance
UPDATE
I want something like this map.

Accepted Answer

Adam Danz
Adam Danz on 8 Feb 2020
Edited: Adam Danz on 8 Feb 2020
The border your plotting (via Greene's borders function), is not the same as the maps produced by the map toolbox. So it's tricky when you mix the two.
Here are some options I came up with that somewhat fit your description. The interpolation only extends as far as the data, not to the latitude/longitude borders.
Produce the base plot
The options that follow will all be applied to this base plot individually.
% Load data
load('EXAMPLE.mat')
% Produce the map
borders('Iran Islamic Republic of') % from FEX
grid on
axis equal % I added this, important to maintain aspect ratio
% ax.GridLineStyle = '-'; % I removed this, it's the default style;
% Modify the X and Y ticks positions
xGrid = 44:.5:65; % I added this because we'll use these values later
yGrid = 25:.5:40; % (same as above)
xticks(xGrid); % Use the variables
yticks(yGrid); % Use the variables
xlabel('Longitude')
ylabel('Latitude')
xtickangle(90)
Show all stations and the rainfall data
hold on
sch = scatter(EXAMPLE.lon, EXAMPLE.lat, 35, EXAMPLE.normalized_rainfall, 'fill','MarkerEdgeColor', 'w');
cbh = colorbar();
ylabel(cbh, ['Rainfall', newline(),'(normalized units)'])
Interpolate the scattered data into a grid
The grid is defined by the axes' grid values. Note that the interpolation will only extend to the extent of the data.
% Interpolate
Vq = griddata(EXAMPLE.lon, EXAMPLE.lat, EXAMPLE.normalized_rainfall, xGrid, yGrid.');
[xq,yq] = meshgrid(xGrid, yGrid);
% plot
hold on
sch = scatter(xq(:), yq(:), 35, Vq(:), 'fill');
cbh = colorbar();
ylabel(cbh, ['Rainfall', newline(),'(normalized units)'])
Create an interpolated 3D surface plot
This uses the exact same interpolated data as above.
% Interpolate
Vq = griddata(EXAMPLE.lon, EXAMPLE.lat, EXAMPLE.normalized_rainfall, xGrid, yGrid.');
[xq,yq] = meshgrid(xGrid, yGrid);
% Plot it
ax = gca();
sfh = surf(xq, yq, Vq, 'EdgeColor', ax.GridColor, 'EdgeAlpha',ax.GridAlpha);
axis tight
cbh = colorbar();
ylabel(cbh, ['Rainfall', newline(),'(normalized units)'])
To have some fun try executing view(3) for a 3D representation; view(2) will return to 2D view.
Also note that you could add the scatter points from the first plot on top of this to see data from single stations.
200208 154012-Figure 1.png
Use imagesc
The advantage of this method over the previous method is that this is 2D instead of 3D and therefore avoids potential problems such as scatter point markers being hidden by 3D peaks even under a 2D view.
The disadvantage of this method over the previous method is that you've got to add in the grid manually since the image consumes the entire plot.
% Interpolate
Vq = griddata(EXAMPLE.lon, EXAMPLE.lat, EXAMPLE.normalized_rainfall, xGrid, yGrid.');
% Plot the gridded image
axh = gca();
hold on
imh = imagesc(xGrid, yGrid, Vq);
axClim = axh.CLim;
axh.Colormap = [1,1,1; axh.Colormap]; % Set NaN values to white
caxis([axClim(1)+axClim(1)*.01, inf]) % This avoids white pixels appearing at the lowest data points
cbh = colorbar();
ylabel(cbh, ['Rainfall', newline(),'(normalized units)'])
% You have to add the border line on top of the image
bh = borders('Iran Islamic Republic of','k-','LineWidth',2);
% Add pseudo grid lines
% note the -0.25 shift to account for imagesc bin centering.
grid off
arrayfun(@(i)xline(xGrid(i)-.25,'-','Color',axh.GridColor, 'Alpha', axh.GridAlpha),1:numel(xGrid));
arrayfun(@(i)yline(yGrid(i)-.25,'-','Color',axh.GridColor, 'Alpha', axh.GridAlpha),1:numel(yGrid));
axis tight
Also note that you could add the scatter points from the first plot on top of this to see data from single stations.
200208 154806-Figure 1.png
  4 Comments
BN
BN on 9 Feb 2020
Thank you, Actually, I chose surf. Yes, I'm done that:
Capture.JPG
It seems that everything is accurate. Thanks for the mention of caxis(), I'll use it.
Thank you again.
Best Regards
Adam Danz
Adam Danz on 9 Feb 2020
It looks like you've layered the gridded scatter points on top as well (hense all of the dots at each grid intersection). That's probably a bit much but at least it shows that the results from all methods agree.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!