Interpolate scattered latitude and longitude (climate) data and plot

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.

[darova]

I tried using griddata

The result:

Answer 1

[Adam Danz]

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.

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.

[BN]

Really Thank you. This solved my problem perfectly.

[Adam Danz]

Glad I could help! Whichever method you use, I would plot the original data from the first scatter plot on top of the interpolated data just to see if the interpolation makes sense and isn't shifted in either direction as a result of binning.

Also you might want to check out caxis() so your color range can extend from -2 to 2 or whatever the normalized limits are.

[BN]

Thank you, Actually, I chose surf. Yes, I'm done that:

It seems that everything is accurate. Thanks for the mention of caxis(), I'll use it.

Thank you again.

Best Regards

[Adam Danz]

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.