Put lat/lon points onto an 80km x 80km grid

Hello,

I have a list of lat/lon points (severe weather reports) that all fall over the continental united states.

What I want to do is basically put these points onto an 80km x 80km grid, and look at whether each grid box has a report in it or not. The number of reports in each box does not matter, just whether it has a report in it or not.

So the end result should be a matrix of 1's and 0's, with each index corresponding to a certain 80km x 80km box over the US.

Thanks

Answer 1

[Kelly Kearney]

The hist3 function will do the main binning for you, once you define your bin edges. A perfect 80km x 80km grid might be a little difficult to set up in terms of lat/lon edges, but assuming you're willing to fudge it a bit, a ~0.719-deg (lat) x ~1.09 (lon) grid will be about the right resolution.

If you have exact coordinates of the 80 km x 80 km grid boxes, inpolygon could help.

[Kyle Chudler]

This is for something scientific so I'd like to be accurate...based on some quick checks, using those set lat/lon grid spacings makes the width of the box vary from 77 to 110 km as laditide decreases. I'm going to need to be more precise than that.

Answer 2

[Chad Greene]

If you need even grid spacing you'll have to work in some projected coordinates. Use projfwd to get your lat,lon data into some projected x,y coordinates. You'll have to choose a projection that minimizes distortion for the area of your data, but for an area the size of the United States every projection will have some distortion.

Check out distortcalc, tissot, and mdistort for different ways to compare map distortions.

After your data are in projected coordinates I'd follow Kelly's advice and use hist3 to determine how many data points are in each grid cell. Using the syntax N = hist3(...) you can then get the logical matrix M you're seeking by M = N>0.

[Kyle Chudler]

This seems like it will work!

However, I'm having some trouble making sense of the projfwd function. Specifically, the proj map projection structure needed. How can I pull this out from a map of the US? I'm currently getting the map just by doing the following:

usamap conus
geoshow(states, 'FaceColor', 'white')