0

I would like to calculate the distance of a portion of a line that is inside a polygon, as well as the portion that is outside of the polygon, based on from/to lat/long points on the line, and the lat/long points that define the polygon boundary.

Use case: A vehicle travels from one to the next. Lat/long points are recorded at each stage of the journey. Lat/long points are also provided defining the state boundary.

How would I determine the distance travelled in each state?

I am really new to this, and I am not using any GIS software. I am hoping to write a function to perform this calculation, but I would really like to understand the math and techniques used to arrive at the answer.

Any references that can be provided would also be really helpful; I don't know all the GIS terminology, so I haven't had too much success with this case on Google so far...

It has been suggested that I use a library to help perform this function, which I am totally open to. However, I still have no idea where to begin with that the methods that should be called, or even what the name of the process is for calculating that. (I'm guessing it's not "distance of a section of a line inside and outside a polygon". Any light that can be shed would be great!

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Andrew Odri
  • 101
  • 3
  • 1
    Since you're wanting to do this programmatically it's not really a duplicate, but it is related to http://gis.stackexchange.com/questions/98518/calculate-line-distances-inside-of-polygon/ and other questions related/linked to that one. The actual question might be answered as part of one or more other questions here already, as I know there are a few that discuss calculating if lines between xy coordinates intersect or if a geometry falls within or intersects another. – Chris W Nov 11 '14 at 22:25
  • Generally speaking I would use a library to do this, I don't know too many people who would do this mathematically. OGR is a good place to start looking, it's open source so all the maths is there for the taking - expect to do a lot of searching though. I would suggest it would be better to use OGR objects to do this - why bark if you've got a dog? – Michael Stimson Nov 11 '14 at 22:25
  • Thanks @ChrisW! I did see that exact question @ChrisW, and even trawled through the related ones, but it didn't really shed too much light on how to actually perform the operation without ArcGIS or even using a library. – Andrew Odri Nov 11 '14 at 22:33
  • Thanks @MichaelMiles-Stimson! I will definitely look deeper into GDAL/OGR Michael; however I still don't know where to begin with the appropriate methods to call for this calculation... I will modify the question shortly to request this as well. Thanks guys! – Andrew Odri Nov 11 '14 at 22:34

1 Answers1

1

There is a good tutorial http://www.gdal.org/ogr_apitut.html

It's written in 3 languages : C, C++ and python (near the bottom)

How to intersect? you get both geometries (OGRgeometry) which you can make from WKT which is easy to write then use OGRGeometry::Intersection to calculate the bit that's inside then cast that to OGRlinestring and call get_length().

When you call get_length() it will return the length in the units that the geometry is in, which given lat/lon points will be in DD (decimal degrees) which is a length, but not particularly useful, to get a length in metres you will need to project the intersection result into a suitable projected coordinate system; projected coordinate systems will have a defined limit of what part of the earth they cover and you can go slightly outside but not too far - the further outside you go the more the shape distorts. A good system is UTM (universe transverse mercator) - note the pictures.

UTM is broken up into zones that are 6 degrees wide (quite a long way, but it is possible when using state sized data to cross more than one), to find out which one to use pick a point on the line (middle would be best) and compare it to the East-West extents of the UTM zones. To project you need to set the known spatial reference (from the GPS) which will most likely be WGS84 (EPSG:4326) using OGRgeometry::assignSpatialReference and then perform the projection using OGRgeometry::transformTo supplying the projected coordinate system.

To create a spatial reference (OGRspatialReference) I would use initFromEPSG as the EPSG codes are easy to calculate... WGS84 Geographic is 4326, UTM (North of the Equator) is 32600 + zone id (UTM 30N is 32630), if you're south of the Equator UTM is 32700 + zone id (UTM 56S is 32756).

If you must have the length in other units (feet, miles, inches, cubits..) then use a scale factor on the determined length in metres... don't go looking for a projected coordinate system with units of feet, miles, inches.. there probably wont be one and there's a possibility that the EPSG code wouldn't be recognized if you did.

Alternately you can use a conic or equal-area projection to cover all of the likely area without too much distortion. Projections exist for Continental United States that covers all the states (except Alaska sometimes and not Hawaii - but you wouldn't be driving out there anyway) which is a catch-all situation. The difference between the methods is likely to be only a few percent and neither takes into account terrain (traversed distance).

Some information on map projections that might help explain what I'm talking about and help decide which to use.

If you are using OGR libraries it's usual to open a layer (covered in the tutorial) and then iterate through the features (OGRfeature) and get their geometries using GetGeometryRef, as for the GPS log it need not exist and may be created in memory using the importFromWKT call.

Michael Stimson
  • 25,566
  • 2
  • 35
  • 74