You have to build the proj.4 string for the CRS from the parameters given:
+proj=lcc +lat_1=44.883333 +lat_2=45.133333 +lat_0=44.791111 +lon_0=-93.383333 +x_0=152400.000000 +y_0=30480.000000 +a=6378418.9409999996 +b=6357033.3098455509 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
Note that x_0 and y_0 have to be in meters, while false Easting and Northing from the definition (and the WKT mkennedy added) are in units of the projection (us-ft).
With that, you can load the data as delimited text, with X
set to the itm
column and Y
to the utm
column. If you have not set QGIS to prompt for the CRS, rightclick on the layer -> Set CRS for layer
to get the right alignment:
I'm not sure why they have named the second column utm, the projection has nothing to do with Universal Tramnsverse Mercator.
This is quick and dirty. Most of this was adapted from the following website:
http://www.movable-type.co.uk/scripts/latlong.html
I haven't tested this too much but it seemed to work after initial testing. It will return a list of lat, long coordinate pairs along a line at a specified interval. I wrote it in python since I don't know php very well.
from math import cos, sin, atan2, sqrt, radians, degrees, asin, modf
def getPathLength(lat1,lng1,lat2,lng2):
'''calculates the distance between two lat, long coordinate pairs'''
R = 6371000 # radius of earth in m
lat1rads = radians(lat1)
lat2rads = radians(lat2)
deltaLat = radians((lat2-lat1))
deltaLng = radians((lng2-lng1))
a = sin(deltaLat/2) * sin(deltaLat/2) + cos(lat1rads) *
cos(lat2rads) * sin(deltaLng/2) * sin(deltaLng/2)
c = 2 * atan2(sqrt(a), sqrt(1-a))
d = R * c
return d
def getDestinationLatLong(lat,lng,azimuth,distance):
'''returns the lat an long of destination point
given the start lat, long, aziuth, and distance'''
R = 6378.1 #Radius of the Earth in km
brng = radians(azimuth) #Bearing is degrees converted to radians.
d = distance/1000 #Distance m converted to km
lat1 = radians(lat) #Current dd lat point converted to radians
lon1 = radians(lng) #Current dd long point converted to radians
lat2 = asin( sin(lat1) * cos(d/R) + cos(lat1)* sin(d/R)* cos(brng))
lon2 = lon1 + atan2(sin(brng) * sin(d/R)* cos(lat1),
cos(d/R)- sin(lat1)* sin(lat2))
#convert back to degrees
lat2 = degrees(lat2)
lon2 = degrees(lon2)
return[lat2, lon2]
def main(interval,azimuth,lat1,lng1,lat2,lng2):
'''returns every coordinate pair inbetween two coordinate
pairs given the desired interval'''
coords = []
d = getPathLength(lat1,lng1,lat2,lng2)
remainder, dist = modf((d / interval))
counter = 1.0
coords.append([lat1,lng1])
for distance in xrange(1,int(dist)):
c = getDestinationLatLong(lat1,lng1,azimuth,counter)
counter += 1.0
coords.append(c)
counter +=1
coords.append([lat2,lng2])
return coords
if __name__ == "__main__":
#point interval in meters
interval = 1
#direction of line in degrees
azimuth = 279.91
#start point
lat1 = 41.97809
lng1 = -91.65573
#end point
lat2 = 41.97834
lng2 = -91.65767
coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
print coords
I hope this helps. sorry if not. :)
Ok, I added the ability for it to calculate the azimuth for you. I tested with the same location you were using and it turned out perfectly.http://www.darrinward.com/lat-long/?id=705140
import math
def getPathLength(lat1,lng1,lat2,lng2):
'''calculates the distance between two lat, long coordinate pairs'''
R = 6371000 # radius of earth in m
lat1rads = math.radians(lat1)
lat2rads = math.radians(lat2)
deltaLat = math.radians((lat2-lat1))
deltaLng = math.radians((lng2-lng1))
a = math.sin(deltaLat/2) * math.sin(deltaLat/2) + math.cos(lat1rads) * math.cos(lat2rads) * math.sin(deltaLng/2) * math.sin(deltaLng/2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
d = R * c
return d
def getDestinationLatLong(lat,lng,azimuth,distance):
'''returns the lat an long of destination point
given the start lat, long, aziuth, and distance'''
R = 6378.1 #Radius of the Earth in km
brng = math.radians(azimuth) #Bearing is degrees converted to radians.
d = distance/1000 #Distance m converted to km
lat1 = math.radians(lat) #Current dd lat point converted to radians
lon1 = math.radians(lng) #Current dd long point converted to radians
lat2 = math.asin(math.sin(lat1) * math.cos(d/R) + math.cos(lat1)* math.sin(d/R)* math.cos(brng))
lon2 = lon1 + math.atan2(math.sin(brng) * math.sin(d/R)* math.cos(lat1), math.cos(d/R)- math.sin(lat1)* math.sin(lat2))
#convert back to degrees
lat2 = math.degrees(lat2)
lon2 = math.degrees(lon2)
return[lat2, lon2]
def calculateBearing(lat1,lng1,lat2,lng2):
'''calculates the azimuth in degrees from start point to end point'''
startLat = math.radians(lat1)
startLong = math.radians(lng1)
endLat = math.radians(lat2)
endLong = math.radians(lng2)
dLong = endLong - startLong
dPhi = math.log(math.tan(endLat/2.0+math.pi/4.0)/math.tan(startLat/2.0+math.pi/4.0))
if abs(dLong) > math.pi:
if dLong > 0.0:
dLong = -(2.0 * math.pi - dLong)
else:
dLong = (2.0 * math.pi + dLong)
bearing = (math.degrees(math.atan2(dLong, dPhi)) + 360.0) % 360.0;
return bearing
def main(interval,azimuth,lat1,lng1,lat2,lng2):
'''returns every coordinate pair inbetween two coordinate
pairs given the desired interval'''
d = getPathLength(lat1,lng1,lat2,lng2)
remainder, dist = math.modf((d / interval))
counter = float(interval)
coords = []
coords.append([lat1,lng1])
for distance in xrange(0,int(dist)):
coord = getDestinationLatLong(lat1,lng1,azimuth,counter)
counter = counter + float(interval)
coords.append(coord)
coords.append([lat2,lng2])
return coords
if __name__ == "__main__":
#point interval in meters
interval = 1.0
#direction of line in degrees
#start point
lat1 = 43.97076
lng1 = 12.72543
#end point
lat2 = 43.969730
lng2 = 12.728294
azimuth = calculateBearing(lat1,lng1,lat2,lng2)
print azimuth
coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
print coords
I think I have fixed the bugs from the original and it should work (including change the interval). Let me know how it turns out.
Best Answer
Load the PDF into QGIS, deselect or remove all the layers your aren't interested in, leaving you with just the contours. You should end up with something like this (using Quick Map Services' OSM base map):
Note that the noise level of these contours is lost. Its an annotation in the PDF that has no geographic reference to the line, so does not load with the other layers from the PDF. So if you need the noise level for each line you need to add a new attribute to each contour line and enter it manually.
Then use
Vector: Geometry Tools: Extract Nodes
to create a Point shapefile of the nodes of the contours. Note these points form the line segments that make the contour lines. If you need more points on the lines then you need to densify the lines first, which is another step.Then with the point layer loaded, it should now look like this:
Now use
Vector: Geometry Tools: Export/Add Geometry Columns
. Decide what coordinate system you want the points in. Don't click the "save to a new shapefile" box, we're going to add them to the current map layer. Hit OK, then when its run hit Close.Now your attribute table should have the X and Y coordinates:
in this case in the original coordinate system.
Saving the shapefile will give you a DBF file that you can read into Excel or other spreadsheet program for further processing and upload to your database...