1. 程式人生 > >【測繪程式設計試題集】 試題10 座標系轉換

【測繪程式設計試題集】 試題10 座標系轉換

座標系轉換

資料

a,6378137.000
1/f,298.3
L0,111

Q71,36.082771,109.191366,5533.025
P91,33.445550,110.154237,4785.906
Q42,38.372964,108.023609,5453.323
Q34,39.305664,111.361612,5800.386
B99,37.264007,108.385066,5557.617
A89,37.371094,112.321633,5282.713
P90,36.035483,111.145753,5550.860
P60,33.334965,109.295619,5496.801
Q89,35.411290,112.303315,5625.113
P72,36.456890,112.47257,4566.878
A05,38.085189,109.182573,4712.183
A46,35.221248,110.24663,5220.175
Q03,36.182447,112.254924,5426.659

問題

在這裡插入圖片描述01
02
03
04
05
06

  1. 橢球類
   package com.te.coorTrans.lib;
   
   import com.sun.org.apache.xml.internal.security.Init;
   import java.util.Map;
   
   /**
    * <p>Title : Ellipsoid </p>
    * <p>Description : 橢球</p>
    *
    * @author huifer
    * @date 2018/11/12
    */
   public class Ellipsoid
{ //長半軸 public double a; //扁率 public double f; //第一偏心率的平方 public double eccSq; //第二偏心率的平方 public double ecc2Sq; public double M0; public Ellipsoid() { a = 6378137.0; f = 0.335281066475e-2; } /*** * 橢球構造 * @param a * @param invf */
public Ellipsoid(double a, double invf) { this.a = a; f = 1 / invf; Init(); } void Init() { double b = a * (1 - f); eccSq = (a * a - b * b) / a / a; ecc2Sq = (a * a - b * b) / b / b; M0 = a * (1 - eccSq); } /*** * 計算W * @param B 緯度(以弧度為單位) * @return */ public double W(double B) { double w = Math.sqrt(1 - eccSq * Math.pow(Math.sin(B), 2)); return w; } /*** * 計算eta * @param B 緯度(以弧度為單位) * @return */ public double Eta(double B) { double eta = Math.sqrt((ecc2Sq * Math.pow(Math.cos(B), 2))); return eta; } /*** * 計算tan * @param B 緯度(以弧度為單位) * @return */ public double Tan(double B) { return Math.tan(B); } /*** * 計算 卯酉圈的曲率半徑N * @param B 緯度(以弧度為單位) * @return */ public double N(double B) { double w = W(B); return a / w; } /*** * 子午圈曲率半徑 M * @param B * @return */ public double M(double B) { double over = a * (1 - eccSq); double w = W(B); return over / (w * w * w); } }
  1. 大地座標 和 笛卡爾座標系轉換
   package com.te.coorTrans.lib;
   
   /**
    * <p>Title : Position </p>
    * <p>Description : 三維座標轉換</p>
    *
    * @author huifer
    * @date 2018/11/12
    */
   public class Position {
   
       private Ellipsoid ellipsoid;
   
       public Position(Ellipsoid ellipsoid) {
           this.ellipsoid = ellipsoid;
       }
   
   
       /***
        * 大地測量轉換為笛卡爾座標
        * @param lat 緯度
        * @param lon 經度
        * @param height 高度
        * @return 笛卡爾座標系 xyz
        */
       public double[] geodetic2Cartesian(double lat, double lon, double height) {
           double X, Y, Z;
   
           double N = ellipsoid.N(lat);
           double slat = Math.sin(lat);
           double clat = Math.cos(lat);
           X = (N + height) * clat * Math.cos(lon);
           Y = (N + height) * clat * Math.sin(lon);
           Z = (N * (1 - ellipsoid.eccSq) + height) * slat;
           return new double[]{X, Y, Z};
       }
   
       public double[] cartesian2Geodetic(double X, double Y, double Z) {
           double lat = 0;
           double lon = 0;
           double height = 0;
           double POSITION_TOLERANCE = 0.0001;
   
           double p, slat, N, htold, latold;
   
           p = Math.sqrt(X * X + Y * Y);
           if (p < POSITION_TOLERANCE) {
               lat = (Z > 0 ? 90.0 : -90.0);
               lon = 0;
               height = Math.abs(Z - ellipsoid.a * Math.sqrt(1 - ellipsoid.eccSq));
           } else {
               lat = Math.atan2(Z, p);
               height = 0;
               for (int i = 0; i < 5; i++) {
                   slat = Math.sin(lat);
                   N = ellipsoid.N(lat);
                   htold = height;
                   height = p / Math.cos(lat) - N;
                   latold = lat;
                   lat = Math.atan2(Z + N * ellipsoid.eccSq * slat, p);
                   if (Math.abs(lat - latold) < 1.0e-9
                           && Math.abs(height - htold) < 1.0e-9 * ellipsoid.a) {
                       break;
                   }
                   lon = Math.atan2(Y, X);
                   {
                       if (lon < 0.0) {
                           lon += 2 * Math.PI;
                       }
                   }
               }
           }
           return new double[]{lat, lon, height};
       }
   
   
   }
  1. 高斯演算法

    1. 正算
    2. 反算
   package com.te.coorTrans.lib;
   
   import java.util.ArrayList;
   import java.util.List;
   
   /**
    * <p>Title : Gauss </p>
    * <p>Description : 高斯演算法</p>
    *
    * @author huifer
    * @date 2018/11/12
    */
   public class Gauss {
   
   
       private Ellipsoid ell;
       // 中央子午線
       private double L0;
       // Y 方向偏移量
       private double Y0 = 500000.0;
   
       /***
        * 構造器
        * @param ellipsoid
        * @param midLon
        */
       public Gauss(Ellipsoid ellipsoid, double midLon) {
           this.ell = ellipsoid;
           L0 = midLon;
       }
   
       /***
        * 高斯正算
        * @param B 大地座標B
        * @param L 大地座標L
        * @return double[]
        */
       public double[] BL2XY(double B, double L) {
           double dl = L - L0;
           double x = 0.0;
           double y = 0.0;
           List<Double> c = coefficient();
           double X = 0;
           X = c.get(0) * B + c.get(1) * Math.sin(2 * B) + c.get(2) * Math.sin(4 * B) + c.get(3) * Math
                   .sin(6 * B)
                   + c.get(4) * Math.sin(8 * B) + c.get(5) * Math.sin(10 * B);
           List<Double> a = coeffA(X, B);
           x = a.get(0) + a.get(2) * dl * dl + a.get(4) * Math.pow(dl, 4) + a.get(6) * Math.pow(dl, 6);
           y = a.get(1) * dl + a.get(3) * Math.pow(dl, 3) + a.get(5) * Math.pow(dl, 5);
           y = y + Y0;
           return new double[]{x, y};
       }
   
   
       /***
        * 高斯反算
        * @param x
        * @param y
        * @return
        */
       public double[] xy2BL(double x, double y) {
           double B = 0.0;
           double L = 0.0;
   
           List<Double> c = coefficient();
           y = y - Y0;
           double Bf1 = x / c.get(0);
           double v = EndPointLat(c, x, Bf1);
           List<Double> b = coeffB(v);
           B = b.get(0) + b.get(2) * y * y + b.get(4) * Math.pow(y, 4) + b.get(6) * Math.pow(y, 6);
           double dl = b.get(1) * y + b.get(3) * Math.pow(y, 3) + b.get(5) * Math.pow(y, 5);
           L = L0 + dl;
           return new double[]{B, L};
       }
   
       /***
        * 計算底點緯度
        * @param coeff 係數
        * @param x x
        * @param Bf1 底店唯獨
        * @return
        */
       private double EndPointLat(List<Double> coeff, double x, double Bf1) {
   
           double X = x, dX = 0;
           double Bf0 = 0;
           do {
               Bf0 = Bf1;
               dX = coeff.get(1) * Math.sin(2 * Bf0) + coeff.get(2) * Math.sin(4 * Bf0)
                       + coeff.get(3) * Math.sin(6 * Bf0) + coeff.get(4) * Math.sin(8 * Bf0)
                       + coeff.get(5) * Math
                       .sin(10 * Bf0);
               Bf1 = (X - dX) / coeff.get(0);
           } while (Math.abs(Bf1 - Bf0) > 0.0000000000001);
           return Bf1;
       }
   
   
       /***
        * 計算係數a陣列
        * @param X 弧度長
        * @param B 維度
        * @return double[] 係數
        */
       private List<Double> coeffA(double X, double B) {
   
           List<Double> aCoeff = new ArrayList<>();
   
           double a0, a1, a2, a3, a4, a5, a6;
           double N, t, eta;
           N = ell.N(B);
           t = ell.Tan(B);
           eta = ell.Eta(B);
           a0 = X;
           a1 = N * Math.cos(B);
           a2 = N * t * Math.pow(Math.cos(B), 2) / 2.0;
           a3 = N * (1 - t * t + eta * eta) * Math.pow(Math.cos(B), 3) / 6.0;
           a4 = N * t * (5 - t * t + 9 * eta * eta + 4 * Math.pow(eta, 4)) * Math.pow(Math.cos(B), 4)
                   / 24.0;
           a5 = N * (5 - 18 * t * t + Math.pow(t, 4) + 14 * eta * eta - 58 * Math.pow(eta, 2) * Math
                   .pow(t, 2)) * Math.pow(Math.cos(B), 5) / 120.0;
           a6 = N * t * (61 - 58 * t * t + Math.pow(t, 4) + 270 * eta * eta - 330 * Math.pow(eta, 2)
                   * Math.pow(t, 2)) * Math.pow(Math.cos(B), 6) / 720.0;
   
           aCoeff.add(a0);
           aCoeff.add(a1);
           aCoeff.add(a2);
           aCoeff.add(a3);
           aCoeff.add(a4);
           aCoeff.add(a5);
           aCoeff.add(a6);
           return aCoeff;
       }
   
   
       /***
        * 計算係數b陣列
        * @param Bf 低點唯獨
        * @return double[]
        */
       private List<Double> coeffB(double Bf) {
           List<Double> bCoeff = new ArrayList<>();
           double Nf = ell.N(Bf);
           double Mf = ell.M(Bf);
           double tf = ell.Tan(Bf);
           double etaf = ell.Eta(Bf);
           double b0, b1, b2, b3, b4, b5, b6;
           b0 = Bf;
           b1 = 1.0 / (Nf * Math.cos(Bf));
           b2 = -tf / (2 * Mf * Nf);
           b3 = -(1 + 2 * tf * tf + etaf * etaf) * b1 / 6.0 / Nf 
            
           

相關推薦

測繪程式設計試題 試題10 座標系轉換

座標系轉換 資料 a,6378137.000 1/f,298.3 L0,111 Q71,36.082771,109.191366,5533.025 P91,33.445550,110.154237,4785.906 Q42,38.372964,108.0236

測繪程式設計試題 試題01

資料 問題 解 #!/usr/bin/env python # -*- coding: utf-8 -*- # @File : Taxi.py # @Author : huifer # @Time : 2018/10/17 20:29 import

測繪程式設計試題 試題02

資料 N 矩陣 10.00 13.50 14.00 13.80 13.90 15.60 13.30 14.50 13.70 14.40 13.50 13.30 15.10 16.40 15.40 14.90 11.30 13.50 17.70 13.30 15.70 14.00

測繪程式設計試題 試題01 計程車軌跡資料計算

資料 問題 解 #!/usr/bin/env python # -*- coding: utf-8 -*- # @File : Taxi.py # @Author : huifer # @Time : 2018/10/17 20:29 impo

測繪程式設計試題 試題02 矩陣卷積計算

資料 N 矩陣 10.00 13.50 14.00 13.80 13.90 15.60 13.30 14.50 13.70 14.40 13.50 13.30 15.10 16.40 15.40 14.90 11.30 13.50 17.70 13.30 15.

測繪程式設計試題 試題03

資料 300,21182.88,-7044.56,14639.48 600,21707.87,-6930.28,13906.68 900,22207.04,-6828.65,13147.66 1200,22679.16,-6738.66,12363.84 150

執行力決定命10《告別做爛好人》

執行力、職場、業務、升職、加薪「音頻原文」http://dwz.cn/6sB3By大家好,我是林琳笨,今天跟大家一起分享《告別做爛好人》這個話題!什麽是爛好人?就是脾氣好,不管對方有沒有道理,是對還是錯,總是和氣生財,不反駁,不拒絕。極品爛好人就是:不拒絕,之後也不做,這種人是最可恨的。 爛好人的結果是什麽呢

網路流24題試題庫(二分圖+最大流)

傳送門     試題庫 I think     點集x,y分別放置試題與型別。源點向x集點連容量為1的邊,x集點向y中其所屬型別連容量為1的邊,y集點向T連容量為所需量的邊,求解最大流若等於總題數

最短路二分圖匹配樹形背包DPDay 10.8

void second eof 最小 span har mes find names T1 最短路 1 #include <cstdio> 2 #include <queue> 3 #include <iostream>

bzoj [JSOI2010]Group 部落劃分 Group二分+並查

group print i++ end cst col sqrt 部落 har 我是zz嗎這麽簡單都寫錯…… 一眼二分,然後判斷的話是枚舉點,然後計算這個點到已有聯通塊的最小距離,如果這個點到一些聯通塊的距離小於當前二分的val,則把這些聯通塊合並起來,這裏用並查集維護,最

概率論與數理統計小結10-1 - 假設檢驗概述

sqrt htm get 依據 事件 http 例如 style 科學 註:終於寫到最激動人心的部分了。假設檢驗應該是統計學中應用最廣泛的數據分析方法,其中像"P值"、"t檢驗"、"F檢驗"這些如雷貫耳的名詞都來自假設檢驗這一部分。我自己剛開進入生物信息學領域,用的最多的就

Chapter4*程式設計總結一*(含原始碼)複製空洞檔案且不把0複製到新檔案

一、寫在前面 不積跬步無以至千里,一點點累積最後達到意想不到的效果。認真對待每一個小細節,一點點改正修訂,往往是問題關鍵所在。 二、coding中遇到的坑 步驟一:建立兩個檔案,一個是空洞檔案,另一個是非空洞檔案,分析比較兩者不同。 2-1 shell中出現亂碼 建立無空

POJ 2771 Guardian of Decency 最大獨立

傳送門:http://poj.org/problem?id=2771 Guardian of Decency Time Limit: 3000MS   Memory Limit: 65536K

軟體包合keil V5 V4 c51 MDK 420-423 953-959 512-526 各版本軟體包下載地址官網

mdk525下載地址: http://az717401.vo.msecnd.net/eval/MDK525.EXE mdk524下載地址: http://az717401.vo.msecnd.net/eval/MDK524.EXE mdk523下載地址: ht

程式設計師都是隱藏的文學家,皮起來你無法招架爆笑合

印象中,蘋果商店裡 APP的更新日誌 通常都是一些簡明扼要 說了等於沒說的描述 但我們終究還是 低估了程式設計師們的才華 他們皮起來你根本想象不到! 我們來看看那些年程式設計師們,那些年笑翻我們的奇葩更新日誌吧 這年頭做一個單身狗 實在是太危險了 連更新個A

模板·並查洛谷 P3367 模板並查

題目:並查集 思路: 複習…… 第一次提交忘寫路徑壓縮T了…… 結論:打過再多遍的模板也要檢查一下啊…… 程式碼: #include<bits/stdc++.h> using namespa

vue-cli 3.x 上線 入坑合

vue-cli 3.x 打包上線 【入坑合集】 問題1:使用vue-cli 3.x開發的專案,開發的時候順利無比,一旦打包上線各種問題就來了。首先是資源裡面報各種各樣的請求錯誤。 解決思路:在src資料夾同級目錄下建立vue.config.js檔案 在vue.config.js

再做食物鏈擴充套件並查

這題很明顯總共有三種關係: 1、x與y同類 2、x吃y 3、y吃x 所以這時候,如果只有一個域就不能表示這3種關係了,比如將(x,y)合併,那麼這時候是x,y哪種關係呢? 所以這時候我們就要將這一個點擴充套件為3個點。 同類域,吃域,被吃域。 具體看程式碼吧: #includ

網路流24題 No.10 餐巾計劃問題 (線性規劃網路優化 最小費用最大流)

【題意】   一個餐廳在相繼的 N 天裡, 每天需用的餐巾數不盡相同。 假設第 i 天需要 ri 塊餐巾(i=1,2,…, N)。 餐廳可以購買新的餐巾,每塊餐巾的費用為 p 分;或者把舊餐巾送到快洗部,洗一塊需 m 天,其費用為 f分;或者送到慢洗部, 洗一塊需 n 天(n>m),其費用為 s<

LCT+並查BZOJ2959[長跑]題解

題目概述 CHNJZ可以在 n 個地方虐場,每次虐場可以踩若干個人。一個地方的人被踩後就不能再踩了(心態已爆炸)。 有 m 個事件:1.地點 x 到地點 y 新建了一條邊。2.地點 x 能踩的人變成了 y 。3.詢問從 x 到 y 最多能踩多少人。 解