[Bf-python] a geoToUTM() function in a bpy module ?

[jerome]

Hello,

I'm currently programming things about city generation for a BGE project 
I have.
open street map is a really valuable input for such need, as you know I 
suppose,
since you can retrieve a lot about city geometry worldwide, and generate 
from it in Blender.

in a recent commit I updated a bit the osm importer to add a better 
projection from lat/lon to blender units.

I think this function, or an equivalent one, should be part of the bpy, 
maybe in mathutils.geometry, or a more suitable location as you wish :

it's really a multi usage function.
this could help to bridge with the osm community,
and with architects too.. for now the ones I know are a bit reluctant to 
Blender but I'm hardly working on it, BGE helps a lot actually.

the tests I'm doing with lxml xml parser are very conclusive to read 
write huge osm or extended osm quickly.
by extended I mean extra tags about height, uvs,utm coords. a kind of 
.bosm format I'm writing.

anyway here's the proposed function, consider it copyleft.
sorry if my proposal does not respect blender guidelines, but I really 
have no time left :s

regards,

Jerome / littleneo


(from math import radians, sin, cos, tan, sqrt)

# given lat and longitude in degrees, returns x and y in UTM (1 KM = 1 
BU ) .
# accuracy : supposed to be centimeter. community feedback needed.
# looks ok so far
# http://fr.wikipedia.org/wiki/Projection_UTM
# http://fr.wikipedia.org/wiki/WGS_84
# http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
# 
http://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/alg0071.pdf
# wiki is your friend (don't ask me about math Im just a writing monkey.)
# jerome.le.chat at free.fr
def geoToUTM(lon, lat) :

     # if abs(lat) > 80 : lat = 80 #wrong coords.

     # UTM zone, longitude origin, then lat lon in radians
     z = int( (lon + 180)  / 6 ) + 1
     lon0 = radians(6*z - 183)
     lat = radians(lat)
     lon = radians(lon)

     # CONSTANTS (see refs.)
     # rayon de la terre à l'équateur
     a = 6378.137
     K0 = 0.9996
     # flattening consts
     f  = 0.0033528106647474805  # 1 / 298.257223563
     e2 = 0.0066943799901413165  # 2*f - f**2
     e4 = 4.481472345240445e-05  # e2**2
     e6 = 3.0000678794349315e-07 # e2**3

     # lat0. 10000 for South, 0 for North
     N0 = 10000 if lat < 0 else 0

     A = (lon - lon0) * cos(lat)
     C = (e2 / (1 - e2)) * cos(lat)**2
     T = tan(lat)**2
     vlat = 1 / sqrt( 1 - e2 * sin(lat)**2 )
     slat = (1-(e2/4)-((3*e4)/64)-((5*e6)/256))*lat - 
(((3*e2)/8)+((3*e4)/32)+((45*e6)/1024))*sin(lat*2) + (((15*e4)/256) + 
((45*e6)/1024) )*sin(lat*4) - ((35*e6)/3072)*sin(lat*6)
     E = 500 + (K0 * a * vlat) * (A + (1-T+C)*((A**3)/6) + (5 - 18 * T + 
T**2) * ((A**5)/120) )
     N = N0 + (K0 * a) * ( slat+vlat*tan(lat)* (A**2/2 + 
(5-T+9*C+4*C**2) * (A**4/24) + (61-58*T+T**2) * A**6/720) )
     return E,N