如何为n kilimeter创建方位方向矩形

Question

假设我有三个四列，比如item、item_纬度、item_经度和方位角

我正在尝试为每个纬度和方位角方向的长度创建一个+/-35度和2公里的矩形

例如。场地x的方位角为45度，因此+35度和-35度至45度为 15度和80度，距离为2公里

我正在尝试获取该区域的所有多边形。我使用了jts.util.GeometricShapeFactory库，但它总是创建一个如下所示的圆。我使用了scala和jts.util.GeometricShapeFactory库

def generatePolygoan(lat:String,long:String,diameterINKM:Double):Option[String] = {
    val latitude = lat.toDouble
    val longitude = long.toDouble
    val diameterInMeters = diameterINKM * 1000
    val shapeFactory = new GeometricShapeFactory()
    //shapeFactory.setNumPoints(64) // adjustable
    shapeFactory.setCentre(new Coordinate(latitude, longitude))
    // Length in meters of 1° of latitude = always 111.32 km
    shapeFactory.setWidth(diameterInMeters / 111320d)
    // Length in meters of 1° of longitude = 40075 km * cos( latitude ) / 360
    shapeFactory.setHeight(diameterInMeters / (40075000 * Math.cos(Math.toRadians(latitude)) / 360))
    val circle = shapeFactory.createCircle()
    Some(circle.toString)
  }

上面的代码创建了围绕lat和long的圆，但是我只需要在方位角+/-35度的方向上创建矩形

请帮帮我，我是几何新手。代码可以在任何编程中，比如JAVA/SCALA/PHP/PYTHON/MYSQL/POSTGRES/MONGODB，甚至可以是数学伪代码

Answer 1

您需要创建一个函数，从另一个点在特定方位角和特定距离上计算一个点（请参见geographic distance and azimuth第1C章），然后调用该函数：

• 从项目点到方位角的循环-距离方位角上的增量
• 在距离处从方位角-增量方位角循环到方位角+增量方位角
• 从距离项目点的方位角+增量循环：

代码：

import math

earthRadius = 6371.0 # earth radius in km

# utility functions working in degrees
def cos(d):
    return math.cos(math.radians(d))

def sin(d):
    return math.sin(math.radians(d))

def arccos(d):
    return math.degrees(math.acos(d))

def arcsin(d):
    return math.degrees(math.asin(d))

def printLatLon(lat,lon,prefix="",sep=" "):
    print(str.format("{}{:10.8}{}{:11.9}",prefix,lat,sep,lon))

# function to compute a point at some distance in some direction
def pointFrom(lat1, lon1, azimuth, distance):
    b = math.degrees(distance / earthRadius)
    a = arccos(cos(b)*cos(90 - lat1) + sin(90 - lat1)*sin(b)*cos(azimuth))
    B = arcsin(sin(b)*sin(azimuth)/sin(a))
    lat2 = 90 - a
    lon2 = B + lon1

# set initial point
p0Lat = 18.3965 # latitude in degrees 
p0Lon = 74.7274 # longitude in degrees 
azimuth = 20.0 # azumuth in degrees 
distance = 2.0 # distance in km (sane unit as earth radius)
deltaAzimuth = 35.0

pList=[[p0Lat，p0Lon]]

# conpute arc from point 0 to azimuth - 35
nbPoints = 5
deltaDistance = distance / nbPoints
print(str.format("delta distance {:.3} km",deltaDistance))
azi = azimuth - deltaAzimuth
for i in range(1,nbPoints+1):
    disti = deltaDistance * i
    #print("disti = "+str(disti))
    piLat,piLon = pointFrom(p0Lat,p0Lon,azi,disti)
    printLatLon(piLat,piLon,prefix=str.format("{:.3} {:.3} : ",azi,disti))
    pList.append([piLat,piLon])

# compute n points point at distance from azimuth - 35 to aximuth + 35
deltaAzPoint  = 5.0
nbPoints = int(2.0  * deltaAzimuth / deltaAzPoint)

print(str.format("delta azimuth {:.3} deg",deltaAzPoint))
az0 = azimuth - deltaAzimuth
for i in range(1,nbPoints+1):
    azi = az0 + i*deltaAzPoint
    #print("az : "+str(azi))
    piLat,piLon =  pointFrom(p0Lat,p0Lon,azi,distance)
printLatLon(piLat,piLon,prefix=str.format("{:.3} {:.3} : ",azi,distance))
    pList.append([piLat,piLon])

# compute points from distance to 0 (p0) for azimuth + 35
nbPoints = 5
print(str.format("delta distance {:.3} km",deltaDistance))
azi = azimuth + deltaAzimuth
for i in range(nbPoints-1,0,-1):
    disti = deltaDistance * i
    piLat,piLon = pointFrom(p0Lat,p0Lon,azi,disti)
    printLatLon(piLat,piLon,prefix=str.format("{:.3} {:.3} : ",azi,disti))
    pList.append([piLat,piLon])

# return to origin
pList.append([p0Lat,p0Lon])

# print result
for latLon in pList:
    printLatLon(latLon[0],latLon[1])