inpolygon function within forloop help

Hi!

I am writing a script that reads in geolocated data, bins the data along a grid, and calculates the mean of the binned data for each grid cell. Essentially, I am wanting to create a raster, but the figure output suggests that I am not binning points in the correct location. If you wouldn't mind, could you take a look at my code and see if I am missing something? I've been working on this for a week now and still haven't made progress.

Thank you in advance!

%Before this code snippet, I load in variables from a csv file x, y, and t.
%x and y are locational data, t is time.

x_bin_edges = (234399:3:235602)';%Vector of x direction bin edges
y_bin_edges = flip(4325799:3:4327002)';%Vector of y direction bin edges
raster_mean = zeros(length(x_bin_edges)-1,length(y_bin_edges)-1);%Allocate a vector space

for i = 1:length(x_bin_edges)-1
    for j = 1:length(y_bin_edges)-1
        in = inpolygon(x,y,[x_bin_edges(i),x_bin_edges(i+1)],[y_bin_edges(j),y_bin_edges(j+1)]);
        if sum(in) < 1%Output NaN's for grid cells that do not have any values.
        raster_mean(i,j) = NaN;
        elseif sum(in) > 0%Output the mean of the data (t) within the grid
        raster_mean(i,j) = mean(t(in));
        end
    end
end

When I plot the data, I get the following:

As you can (hopefully) see, the red +'s are the plotted locations of the original data, whereas the colored grid cells is the plotted raster_mean calculated in the for loop. Not only do the data not line up, but the plot seems to suggest missing information.

[darova]

I don't understand the question. Can you please explain more? Can you make a sketch or something?

[Randall Bonnell]

Hi darova!

Thank you for the comment. I'm sorry for the lack of clarity. Perhaps it would be helpful for me to better explain the figure?

Here's the code I used to create the figure:

figure(4)%Plot the raster_mean results on top of original data
clf
hold on
mapshow(x,y,'DisplayType','point')
mapshow(raster_mean,R,'DisplayType','Surface')

Here, the R is a MapCellsReference Structure created by the readgeoraster function. It's bounding limits and resolution are what determined the bin edges for the forloop I used. I am trying to create a matrix using the forloop that will plot values over the original data shown as red +'s in the figure above. As you can see in the figure, the values do not plot over the original data.

Does that make more sense?

[darova]

Can you show the picture you want to get?

Something like this

[Matt J]

@Randall BonnellBut raster_mean consists of average t-values, not x,y locations. How can you therefore form a scatter plot of them?

[Randall Bonnell]

Hi @Matt J,

Thank you for your comment! Prior to these lines of code, I read in a tiff:

[M,R] = readgeoraster('data.tif','OutputType','double');
northing_lim = [4325800 4327000]; 
easting_lim = [234400 235600];
[M,R] = mapcrop(M,R,easting_lim,northing_lim);

I then use the R, which is a map cells reference structure to map raster_mean.

...I am also realizing that the limits in the northing_lim and easting_lim do not match the limits set by x_bin_edges and y_bin_edges...That may actually be my issue.

[Randall Bonnell]

Hi @darova,

Thank you for your help on this question! I apologize for my non-specificity.

Answer 1

[Matt J]

I don't see anything wrong with your code, but this might be simpler:

m=numel(x_bin_edges)-1; n=numel(y_bin_edges)-1;

[~,~,~,binX,binY] = histcounts2(x,y,x_bin_edges,y_bin_edges);

raster_mean=accumarray([binX,binY],t,[m,n],@mean);

[Randall Bonnell]

Hi again @Matt J,

Thank you so much for this suggestion. I'll attempt to use it and see if it improves the run-time.

[Randall Bonnell]

Hi @Matt J,

That did the trick. Thank you so much for your help, it seemed like the inpolygon function wasn't binning data along the edges of the bins.