How to automate selecting specific xyz map data from an existing xyz data set?

I have a bathymetry map that provides a 1 meter resolution. However, I don't want all of the map data. I want a 1024X1024 meter section of the map so that I can run comparrisons. The xyz map data I have is lat/long/alt and I'd like to automate the process and be able to populate a 1024X1024 matrix with z data.

I have attached a link to the xyz file since it is too large and a .jpg file so that you can see what the map looks like. I've also drawn a blue box around the desired location.

http://seafloor.otterlabs.org/SFMLwebDATA_mb.htm#MONT

This is a painful process to go cell by cell so I appriciate all help.

[Walter Roberson]

The section you are extracting: is it a cuboid aligned with linear coordinates? A cuboid aligned with lat/long coordinates? A cuboid at an angle to linear or lat/long coordinates? An arbitrary shape?

[BlkHoleSun]

First, thanks for responding!

It is aligned with lat/long coordinates,

[Walter Roberson]

This would be a little easier if you used the Shaded Relief Images that are just a couple of links higher up than what you circled. Those are in geotiff form, and the Mapping Toolbox has tools to extract sections from geotiff.

[BlkHoleSun]

I am open to an easier method. i just used the [A R]= geotiffread(filename) function to extract the data. Here is R:

[Walter Roberson]

Bleh. The tif file is a map reference rather than a geo reference, so it is projected coordinates into some system. Unfortunately the information present there does not appear to talk about what the coordinates mean.

When I load the GeoTiffs directory and look inside the .xml files, I see that it is Transverse Mercator centered around long -123.0 lat 0. There is transform information there but I do not know how to understand it yet.

It looks like it might be Geodetic Reference System 80 (GRS80) ... ah, but looking at other items it looks like WGS84 is also there.

Hmmm, and there that information present using geotiffinfo() projection parameters sub-structure.... But I still do not understand yet how to convert to lat/long. Ah, I made some progress at http://www.rcn.montana.edu/resources/converter.aspx treating the readings as UTM in zone 10N...

(I have done very little map work, so I'm learning as I go along.)

[BlkHoleSun]

I've made some good progress so far. I've been able to get the lat/lon values and create a box around what I want. I was attempting to extract the data by finding the values that are > or < the lines I made around the box. That doesn't seem to be working since the find function doesn't like non-constant threshold values. Any thoughts or advice on extracting the desired data?

[Walter Roberson]

A question would be how accurate the boundaries have to be. The map spacing is in physical meters, but imposing lat/long boundaries implies boundaries that are rectilear in physical meters. Your source about 7.3 km by about 6.9 km, which I figure is large enough that the curvature of the Earth would several few meters off of the upper left and upper right if the lat/long boundaries need to be strict.

I would tend to suggest that you just translate the lat/long boundaries at the southern edge and use those to select a rectangular portion of the data instead of worry about the curvature.

Answer 1

[BlkHoleSun]

Walter and KSSV,

Sorry it's been a while. I figured it all out. I first determined the slope of the box I wanted and the top left point (could have pick any point). I then used algebra to compute the other corners of the box. I then did a linspace to make 2 sets of parallel lines with 256 points. I then created a for loop to connect the two lines using linspace again for 256 points inbetween the two lines. This resulted in a 256X256. Next, I had to get z data over. I used the reshape function to get my 256X256 into a single column, then used the griddata function for the z values and used reshape to get it back into a 256X256. Now I'm on to calculating Boltsman and Shannon entropy measure for comparrison! Thanks for all your help!

Answer 2

[KSSV]

Read about inpolygon. If you have the vertices of closed region/ polygon; ou can get the points lying inside the region.

[Walter Roberson]

So far the boundaries are expected to be in latitude and longitude, but the data is in Transverse Mercator Projection (basically an unrolled local cylindar.) Effectively the data is in ground distance and lat and long boundaries are curves relative to that. You would not be able to use a rectangular polygon for inpolygon purposes: you would have to approximate the lat/long boundaries by a series of straight lines. A bit of a nuisance at the very least.

[BlkHoleSun]

wow! this was a super helpful function. I was able to quickly identify all points inside of the box. I only have one more question.

How do I get the altitude values that correspond to the x and y points inside of the square. Those values are in a seperate matrix.

Thanks!

[KSSV]

You need to do interpolation. Read about griddata and/ or scatteredInterpolant.

Or alternatively, you have the indices of the points in hand...you can striaght away pick the respective z values.

[BlkHoleSun]

Walter,

I used a series of lines using some old algebra and approximated the four poins of the square. It was a huge pain.

KSSV,

I'll look into griddata in the morning and let you both know how the results fair.

thanks again for all your help so far...this is a pretty big project!

[KSSV]

YOu are welcome..