How to grid data with coordinates to create a spatial plot using geoshow

I have 3 vectors: data, lat, lon which I am trying to plot spatially for the continental US. Is there a function where I can organize my lat and lon vectors in the appropriate gridded format which geoshow will plot properly while ensuring the data vector is organized in the same fashion so the data points remain with their respective coordinates?

Answer 1

[Chad Greene]

Is it possible that your data vectors are in fact regular, but simply not in gridded format? To check, try

scatterm(lat,lon,20,data)

If the datapoints are regularly spaced, you should be able to use my xyz2grid function like this:

[LON,LAT,DATA] = xyz2grid(lon,lat,data);

where lowercase are your data vectors and upper case are gridded. Then you can do

pcolorm(LAT,LON,DATA)

[Ronnie Abolafia-Rosenzweig]

unfortunately the data vectors are not regular; the scatterm function produced an error when running the program

[Chad Greene]

What was the error message?

[Ronnie]

Error using gcm (line 25) Not a map axes.

Error in scatterm (line 48) gcm(ax);

[Chad Greene]

That means a map hasn't been initialized. Try

worldmap([min(lat) max(lat)],[min(lon) max(lon)]) 
scatterm(lat,lon,20,data)

[Ronnie]

Thank you! This is working very well. I see the 20 denotes the area of each marker; do you know what units this is in? (i.e. if each point is symbolic of a 20 km x 20 km grid cell should I use 20?)

[Chad Greene]

The 20 isn't really units of real-world space. Rather, it's more like the number of pixels it takes up on your screen.

The scatterm suggestion was just to check whether your data correspond to a regular grid. Do the circles produced by scatterm make a gridded pattern? If yes, I recommend

% Grid the data: 
[LON,LAT,DATA] = xyz2grid(lon,lat,data); 

% Initialize a map: 
worldmap([min(lat) max(lat)],[min(lon) max(lon)]) 

% plot the gridded data: 
pcolorm(LAT,LON) 
shading interp

[Ronnie]

I see what you are saying. Unfortunately my arrays are too large for the xyz2grid function to work. I am getting the following error message:

Error using accumarray
Requested 3857x5166498 (148.5GB) array exceeds maximum array size preference. Creation of
arrays greater than this limit may take a long time and cause MATLAB to become
unresponsive. See array size limit or preference panel for more information.

Error in xyz2grid (line 92)
Z = accumarray([yi xi],z(:),[],[],NaN);

Error in Monthly_Average_Function (line 85)
[X,Y,Z] = xyz2grid(double(Z),double(lon),double(lat));

I went into my preferences to uncheck the box that limits the array size limit to RAM; however once I did this MATLAB just shuts down. Do you have any suggestions on this issue?

[Chad Greene]

Oh wow, that's a lot of data points. I haven't downloaded the data you shared, but from KSSV's comment it looks like he found a lot of data points you can get rid of. If it's true, perhaps you can follow his suggestion

%%remove -9999
data(lon==-9999) = [] ;
lon(lon==-9999) = [] ;
lat(lat==-9999) = [] ;

And then use xyz2grid.

[Ronnie]

Even after removing these values I still get the same error message. Is there an advantage to using pcolor vs. scatterm ? I just want to ensure the spatial plot I am creating is an accurate representation

[Walter Roberson]

pcolor and scatterm are quite different purposes.

pcolor does a surface plot of a data grid -- it is surf() followed by view(2). pcolor of an N x M data grid produces (N-1) x (M-1) faces with the colors being interpolated by adjacent data points.

scatterm does a 2D scatter plot on a mapping axes, putting a marker at each location specified.

In answer to your earlier question, the area of the marker, such as 20, is in "points squared" where 1 point is 1/72 of an inch. In the case of circular markers, for size S, you would take R = sqrt(S/pi) to get the radius of the circle, in points. This size is used for the marker no matter how far in or out you zoom, so it does not reflect any particular data units.

[Ronnie]

Thank you, Chad and Walter. I greatly appreciate your insights and I will move forward and utilize this knowledge.

Best, -Ronnie

[Ronnie]

Update: I got it to work by trimming down the data even further (getting rid of NaN values) and then xyz2grid worked; and I am able to show the data using geoshow. Thank you so much for your help!

[Chad Greene]

Great, I'm glad you got it working!

Answer 2

[KSSV]

Let data be your nx3 array which has lon, lat and data in the first, second and column respectively.

    % Get longitude and latitude vectors 
    x = unique(data(:,1)) ;
    y = unique(data(:,2)) ;
    % dimensions of the data
    nx = length(x) ; 
    ny = length(y) ;
    % Frame matrix of grid 
    D = reshape(data(:,3),[ny,nx]) ;
    % flip  matrix to adjust for plot
    H = flipud(H) ;
    % Transpose the matrix 
    H = H' ;  % Check if is required
    surf(x,y,H) ;

[Ronnie Abolafia-Rosenzweig]

In this solution, using the unique function to create latitude and longitude vectors will erase any duplicates of latitude and then any duplicates of longitude. It is okay to have multiple latitude readings, as long as they are paired with different respective longitude readings. I do not want to erase data points from the vectors. Both vecotrs (lat and lon) have already been updated to ensure that there are no duplicate coordinates.

I am seeking for a way to transform vectors in matrices which can be read by the function geoshow.

[KSSV]

The above code works..if your data is a structured grid and in (x,y,z) format......if it is unstructured grid, you need to follow the below:

x0 = min(lon) ; x1 = max(lon) ; nx = 100 ;
y0 = min(lat) ; y1 = max(lat) ; ny = 100 ;
x = linspace(x0,x1,nx) ;
y = linspace(y0,y1,ny) ;
[X,Y] = meshgrid(x,y) ;
Z = griddata(lon,lat,data,X,Y) 
surf(X,Y,Z)

[Ronnie]

My data is not a structured grid. And I cannot use meshgrid for my lat and lon vectors unfortunately because the vector sizes I am using gives me an error saying it is too large to use in meshgrid. I appreciate your help!

[Ronnie]

I am also looking to do a 2d plot, thus the surf function will not work, which is why I am trying to use geoshow.

[KSSV]

Attach your data...it shall work...

[Ronnie Abolafia-Rosenzweig]

I appreciate your help. I will send the data tomorrow once I am back at the computer!

[Ronnie]

SMAP_March was the 3-columned matrix with data in column 1 and the data's corresponding lat and lon in columns 2 and 3 respectively.

data=SMAP_March(:,1); lat=SMAP_March(:,2); lon=SMAP_March(:,3);

%Saves a variable as the script file (will replace rest of script file) matlab.io.saveVariablesToScript('DATA.m','data') matlab.io.saveVariablesToScript('LAT.m','lat') matlab.io.saveVariablesToScript('LON.m','lon')

I used the following code to save the 3 vectors I am sending to you (attached).

[KSSV]

Week end...I will get back to you on Monday ..:)

[Ronnie Abolafia-Rosenzweig]

Enjoy the weekend, I will not have access to my computer until Monday as well. :)

[KSSV]

Where is .mat file?

[Ronnie Abolafia-Rosenzweig]

The site will not allow me to send any files greater than 5 MB, which the lat and data files are. The LON.mat file was bellow the size limit and is attached. Can you please send me an email at [e-mail] and I will respond with the .mat files. I apologize for the inconvenience

[KSSV]

You may attach the file in your google drive and past the link here......

[Ronnie Abolafia-Rosenzweig]

Here is the link to the data: https://drive.google.com/open?id=0B5NUi1TZ3RqcZ2phOFE3TVhVQjA

It should allow you full access. Thank you

[KSSV]

This works....

lon = load('LON.mat') ;
lon = lon.lon ;

lat = load('LAT.mat') ;
lat = lat.lat ;

data = load('DATA.mat') ;
data = data.data ;

%%remove -9999
data(lon==-9999) = [] ;
lon(lon==-9999) = [] ;
lat(lat==-9999) = [] ;

scatter(lon,lat,data,data)

[Ronnie]

Thank you for the response! Much appreciation.

-Ronnie

[TAPAS]

The code xyz2grid is not working it's showing mistake in line 31 in xyz read and line 72 in in xyz2grid