# How do you produce a output table containing distances between mineral grains at various azimuthal directions?

I am using ArcGIS 10, the newest version and I am a senior at Western Carolina University and working on my senior thesis and the data set in which I have are images of mineral thin sections all separated into different groups of grains, magnetite, quartz, etc.. either being a vector layer or a raster layer. I am wanting to calculate the distance between each grain or the grain spacing, boundary to boundary, NOT center to center and also have those distances given every 10 degrees on a 0 to 360 degree directional scale. Is there any code out there that does this very exact method? Or is there another program, R or Fragstats, that would do a better job, ArcGIS 10 preferably.

An example of the grains is shown here:

Thanks,

-Michael Stillwagon

-
You should add that you sent this same question into the R-sig-Geo group, stat.ethz.ch/pipermail/r-sig-geo/2011-February/010916.html – Andy W Feb 14 '11 at 13:57
Is the grain dataset small enough to fit entirely into memory? – Kirk Kuykendall Feb 14 '11 at 15:14
@user1943 Could you please explain what it means to have "distances given every 10 degrees"? One guess is that you want distances binned by relative orientation between grains, but this is ill-defined because the bearing from one grain to another is not well defined. Another guess (based on comments to some existing replies) is you want to restrict the calculation to 1D transects passing through grain centroids, but then do you want actual 2D grain-to-grain distances or do you want the distances along the transects? – whuber Feb 14 '11 at 18:57
@whuber I want to calculate the distance from one grain boundary to the next grain boundary on a angled transverse. From the origin of each grain I want to project an angle every 10 degrees from 0 degrees to 360 degrees this feature will be a line I assume. I want to then measure the distance along each of those "angled transverses" the distance to the nearest grain boundary. This output would have a table containing distances to the nearest grain boundary from a specific angle given. – user1943 Feb 14 '11 at 20:34

To answer one small part of your question:

Using, CAD or similar I would create a template of a circle with radii drawn at the angle you specify. The I would look at where the lines intersect the polygon (see ArcGIS help on intersect).

Then you can use these points for your distance calculations.

As an aside, if you could include a sketch diagram of what you are proposing it would help. For example, it's not clear how many distance measurements you want, or how many 'grains' you have.

-
I do not know how to import an image on here so I have an image with 40 irregular shaped grains having sharp edges and the grains vary in size from 1 micron to .5 micron and the distances between them all vary some being so close that they seem to touch and some being a grain width apart from each other. I want to calculate the distance to the nearest grain on a transverse of a specific angle, so on the 10 degree transverse from the center of the grain I want to know the distance to the nearest grain and at 20 degrees I want to know the distance to the nearest grain and so on. – user1943 Feb 14 '11 at 2:26
Could you post it on imgur.com and I will place in your post? – djq Feb 14 '11 at 2:32
Here is the thin section image in which I am working with: i.imgur.com/o3Zyj.jpg – user1943 Feb 14 '11 at 3:00
just added it. I think it would help if you had some annotations showing your criteria (in particular scale - I did not realize it was this small on first reading it). I know you describe it in detail, but graphics help. And you'll probably need such a diagram for your thesis anyway. – djq Feb 14 '11 at 3:12

I was going to say something similar to celenius, I imagine the workflow being

``````1. get centroid of polygon
2. create polygon ray in direction
3. clip mineral section polygons with ray polygon
4. get distance matrix with those clipped polygons (boundary to boundary)
``````

This could probably be simplified if you knew how to code the nearest neighbor distance from boundary to boundary, so this is assuming you rather just use a package that already has this functionality. I'm sure this can be done in most language environments, but some potential R libraries to do this are listed below.

Various R packages offer the functionalities of 1, 2, and 3. Upon my search it appears the PBSmapping and the sp packages are good options to perform these three functions (with the same objects). I couldn't find any documentation that suggests either package could do 4, but it appears the cshapes package may be able to accomplish this (it is not clear to me if this function will only work on their data or on other data coerced into the same format). The cshapes article I linked to just says it is a wrapper for JTS, so it may be simpler just to do this directly in that environment than worry about R.

If you do code a solution I ask you post it back here on the forum. I'm sure it will be of interest to the general community. I also don't doubt other software packages offer simpler solutions than R, I only list R as I am more familiar with it than other languages.

-

Hallo

Do I understand you right :

1) You want for each grain 36 distances to the surrounding grains. One from each sector of 10 degrees, is that right?

2) Those sectors you are defining from the center of the grain (from a comment to celenius)

3) You want the boundary to boundary distance

If I am right so far I think thee is one thing to comment.

Since you are defining the sectors from the center of the grain and want the boundary to boundary distance it will be kind of messy. If we start to do this calculation from grain1 we will first find the center of it, then define the sector and make a distance calculation from the boundary of grain1 to the boundary of the other grains that falls inside the sector. That means that you might get the shortest distance from a point outside your sector on the boundary of grain1 to a point on the boundary inside the sector on the other grains.

This is doable in PostGIS for instance.

But it would make more sense to me to pick on of the following instead (of cpource depending what you are trying to get out of this)

a) take the distance from the center of grain1 to the boundary inside the sector of the other grains
b) Calculate the distance between the segment of the boundary falling inside each sector both for grain1 and the other grains (means you have to extract also the sector part of the boundary for grain1)
c) Find the "shortest line" between grain1 and the others that satisfy being in each sector

The last one is the most tricky I don't know if it is possible in any GIS-software of today. But in PostGIS it would be quite simple to implement I think. It would be a distance-function with azimuth-restrictions. It would beside the distance also give the coordinates from where the restricted distance is found. That is the beauty with open source. Everything is possible :-) It would need maybe a day or two to get a working patch and if the community likes it, it would be official by summer in PostGIS 2.0.

If you want to stay with your approach as I have understood it above I can post some SQL-code for PostGIS. (I just had to try it when I saw the question :-) )

Approach a) and b) is also doable quite simple.

And as mentioned also approach c) of course also possible

Regards

/Nicklas

Update:

Ok, here comes some code. It is a little bit messy that's why I didn't post it at first. There is some parts that I see could be optimized, but it seems to work anyway.

What it does is: 1) finding the centroid of the first geometry 2) making a buffer around that point 3) finding 36 points around the boundary of the circle from the buffer 4) builds slices from those points to the center point 5) Makes an intersection between the slices and the surrounding grains 6) calculates the distance to those intersections.

That's about it, here comes the SQL:

``````CREATE TABLE result_table AS
SELECT DISTINCT ON (id1, n) id1, dist, id2, n, ST_Shortestline(ST_Intersection(slice, unsliced_geom),  orig_geom) AS shortestline FROM
(
SELECT ST_Distance(ST_Intersection(slices.geom, t1.geom),  orig_geom) AS dist, slices.geom AS slice, t1.geom as unsliced_geom, slices.id as id1, t1.id as id2, slices.n,orig_geom FROM
(
SELECT ST_MakePolygon(ST_LineFromMultiPoint(ST_Collect(array[the_p, p1, p2, the_p]))) AS geom,the_p, id,orig_geom,  n FROM
(
SELECT ST_Line_Interpolate_Point(the_geom, ((n-1)/36.0)::float) AS p1,ST_Line_Interpolate_Point(the_geom, (n/36.0)::float) AS p2,the_p, orig_geom, id, n FROM
(
SELECT ST_Boundary(ST_Buffer(the_p, 1000)) AS the_geom, the_p, id, orig_geom FROM
(
SELECT ST_Centroid(geom) AS the_p, id, geom AS orig_geom FROM the_table
) p
) lw,
generate_series(1,36) n
) points
) slices
INNER JOIN
the_table as t1
ON ST_Intersects(slices.geom, t1.geom) AND slices.id!=t1.id
) t

where id1=1
order by id1,n , dist ;
``````

The WHERE clause makes it only compute against row with id 1. You can just remove that row if you want it to calculate the same thing from each grain.

the resulting shortest line from a small test case:

It can be modified to measure from the sliced part of the "orig_geom" too.

Regards

/Nicklas

-
Nicklas, Thank you for the advice, if you wouldn't mind please post some SQL-codes for PostGIS I would like to try this method and see the results. Based on the information you provided this is certainly doable. Thanks again. – user1943 Feb 14 '11 at 13:50