{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"lunana  \nlast update 2022 05 10  \nゆっくりしていってね","metadata":{"id":"g7Op95I1I1Rv"}},{"cell_type":"markdown","source":"# Data list  \n* [**train**](#train)  \n* [**test**](#test)\n* [**metadata**](#metadata)\n* [**sample_submission.csv**](#sample_submission.csv)","metadata":{"id":"1SN8EgzMUZDo"}},{"cell_type":"markdown","source":"**霊夢:今日はgoogleのコンペだね。  \n魔理沙:まずは概略を見るぞ。**","metadata":{"id":"NtXMQxH4I93e"}},{"cell_type":"markdown","source":"大会の目標  \nこのコンペティションの目標は、スマートフォンの位置をデシメートルまたはセンチメートルの解像度まで計算することです。これにより、HOVレーンETA推定などのレーンレベルの精度を必要とするサービスが可能になります。ホストが収集したデータセットを使用して、オープンスカイおよびライトアーバンロードで収集されたAndroidスマートフォンからの生の位置測定に基づいてモデルを開発します。\n\nあなたの仕事は、より正確な位置を生成するのに役立ち、より細かい人間の行動の地理空間情報と改善された粒度のモバイルインターネットとの間の接続を橋渡しします。その結果、より正確なデータに基づいて新しいナビゲーション方法を構築できます。\n\nコンテクスト  \n高速道路の出口の前で車線変更を逃したことがありますか？他の車線ではなく、相乗り車線の到着予定時刻（ETA）を知りたいですか？これらおよびその他の便利な機能には、正確なスマートフォンの測位サービスが必要です。機械学習モデルは、グローバルナビゲーション衛星システム（GNSS）データの精度を向上させることができます。より洗練されたデータがあれば、何十億ものAndroid携帯ユーザーがより微調整されたポジショニング体験をすることができます。\n\nGNSSチップセットは、スマートフォンの位置を計算するために使用できる生の測定値を提供します。現在の携帯電話は、3〜5メートルの測位精度しか提供していません。高度なユースケースの場合、結果は十分に細かくなく、信頼性もありません。都市の障害物は、GPSの精度に対する最大の障壁となります。このチャレンジのデータには、オープンスカイとライトアーバンロードで収集されたトレースのみが含まれます。これらの高速道路とメインストリートは最も広く使用されている道路であり、スマートフォンのポジショニングの限界をテストします。\n\nGoogleのAndroidGPSチームは、2021年にスマートフォンデシメートルチャレンジを主催しました。3人の受賞者の作品は、ION GNSS+2021カンファレンスで発表されました。今年は、Institute of Navigationの共催で、このコンテストはスマートフォンのGNSS測位精度の高度な研究を求め続け、人々が周囲の世界をよりよくナビゲートできるように支援します。昨年の進捗状況に基づいて構築するために、データには2021年の競争からのトレースも含まれています。\n\n将来の競技会には、衛星信号の障害物がある深部都市部など、より過酷な環境で収集された痕跡が含まれる可能性があります。この競争でのあなたの努力は、このより難しいデータがどのように解釈されるかに影響を与える可能性があります。デシメートルレベルの位置精度により、モバイルユーザーは、より優れた車線レベルのナビゲーション、ARウォーク/ドライブ、電話による正確な農業、および交通安全の問題の場所のより高い特異性を得ることができます。また、よりパーソナライズされた微調整されたナビゲーションエクスペリエンスが可能になります。  \n\nGoal of the Competition  \nThe goal of this competition is to compute smartphones location down to the decimeter or even centimeter resolution which could enable services that require lane-level accuracy such as HOV lane ETA estimation. You'll develop a model based on raw location measurements from Android smartphones collected in opensky and light urban roads using datasets collected by the host.\n\n\n\nYour work will help produce more accurate positions, bridging the connection between the geospatial information of finer human behavior and mobile internet with improved granularity. As a result, new navigation methods could be built upon the more precise data.\n\nContext  \nHave you ever missed the lane change before a highway exit? Do you want to know the estimated time of arrival (ETA) of a carpool lane rather than other lanes? These and other useful features require precise smartphone positioning services. Machine learning models can improve the accuracy of Global Navigation Satellite System (GNSS) data. With more refined data, billions of Android phone users could have a more fine-tuned positioning experience.\n\nGNSS chipsets provide raw measurements, which can be used to compute the smartphone’s position. Current mobile phones only offer 3-5 meters of positioning accuracy. For advanced use cases, the results are not fine enough nor reliable. Urban obstructions create the largest barriers to GPS accuracy. The data in this challenge includes only traces collected on opensky and light urban roads. These highways and main streets are the most widely used roads and will test the limits of smartphone positioning.\n\n\n\nThe Android GPS team in Google hosted the Smartphone Decimeter Challenge in 2021. Works by the three winners were presented at the ION GNSS+ 2021 Conference. This year, co-sponsored by the Institute of Navigation, this competition continues to seek advanced research in smartphone GNSS positioning accuracy and help people better navigate the world around them. In order to build upon last year’s progress, the data also includes traces from the 2021 competition.\n\nFuture competitions could include traces collected in harsher environments, such as deep urban areas with obstacles to satellite signals. Your efforts in this competition could impact how this more difficult data is interpreted. With decimeter level position accuracy, mobile users could gain better lane-level navigation, AR walk/drive, precise agriculture via phones, and greater specificity in the location of road safety issues. It will also enable a more personalized fine tuned navigation experience.","metadata":{"id":"BjDh8x4vKbC9"}},{"cell_type":"markdown","source":"# train ","metadata":{}},{"cell_type":"markdown","source":"**霊夢:次はtrainデータを見てみよう**  \n\n**Reimu: Let's look at the train data next**","metadata":{"id":"eVb1uh75KejL"}},{"cell_type":"markdown","source":"[train / test] / [drive_id] / [phone_name] /device_gnss.csv-各行には、生のGNSS測定値、派生値、およびベースライン推定位置が含まれます。このベースラインは、correctedPrMと衛星位置を使用して、標準の加重を使用して計算されました。最小二乗（WLS）ソルバー。電話の位置（x、y、z）、クロックバイアス（t）、および各エポックの状態としての一意の信号タイプごとのisrbMを使用します。一部の生の測定フィールドは、非推奨であるか、元のgnss_log.txtに入力されていないため、このファイルに含まれていません。\n\n* MessageType-文の接頭辞「Raw」。\n \n* utcTimeMillis-GnssClockから変換されたUTCエポック（1970/1/1）からのミリ秒。\n\n* TimeNanos-GNSS受信機の内部ハードウェアクロック値（ナノ秒単位）。\n\n* LeapSecond-時計の時刻に関連付けられたうるう秒。\n\n* FullBiasNanos-GPS受信機内のハードウェアクロック（getTimeNanos（））と、1980年1月6日の0000Z以降の実際のGPS時間とのナノ秒単位の差。\n \n* BiasNanos-クロックのサブナノ秒バイアス。\n\n* BiasUncertaintyNanos-クロックのバイアスの不確実性（1シグマ）（ナノ秒単位）。\n \n* DriftNanosPerSecond-時計のドリフト（ナノ秒/秒）。\n \n* DriftUncertaintyNanosPerSecond-時計のドリフトの不確実性（1シグマ）（ナノ秒/秒）。\n \n* HardwareClockDiscontinuityCount-ハードウェアクロックの不連続性の数。\n \n* Svid-衛星ID。\n \n* TimeOffsetNanos-測定が行われた時間オフセット（ナノ秒単位）。\n \n* State-衛星の同期状態を示す整数。整数の各ビットは、測定の特定の状態情報に帰属します。ビットと状態の間のマッピングについては、 metadata/raw_state_bit_map.jsonファイルを参照してください。\n \n* ReceivedSvTimeNanos-受信したGNSS衛星時刻（測定時）（ナノ秒単位）。\n \n* ReceivedSvTimeUncertaintyNanos-受信したGNSS時間の誤差推定（1シグマ）（ナノ秒単位）。\n \n* Cn0DbHz-dB-Hz単位の搬送波対雑音比。\n \n* PseudorangeRateMetersPerSecond-タイムスタンプでの疑似距離レート（m / s）。\n \n* PseudorangeRateUncertaintyMetersPerSecond-疑似範囲のレートの不確かさ（1-シグマ）（m / s）。\n \n* AccumulatedDeltaRangeState-これは、「累積デルタ範囲」測定の状態を示します。整数の各ビットは、測定の状態に帰属します。ビットと状態の間のマッピングについては、 metadata/accumulated_delta_range_state_bit_map.jsonファイルを参照してください。\n \n* AccumulatedDeltaRangeMeters-最後のチャネルリセット以降の累積デルタ範囲（メートル単位）。\n \n* AccumulatedDeltaRangeUncertaintyMeters-累積デルタ範囲の不確実性（1シグマ）（メートル単位）。\n \n* CarrierFrequencyHz-追跡された信号のキャリア周波数。\n \n* MultipathIndicator-イベントの「マルチパス」状態を示す値。\n \n* ConstellationType-GNSSコンステレーションタイプ。人間が読める値へのマッピングは、metadata/constellation_type_mapping.csvファイルで提供されます。\n \n* CodeType-GNSS測定のコードタイプ。最近のログでのみ利用可能です。\n \n* ChipsetElapsedRealtimeNanos-システムの起動からこのクロックの経過リアルタイム（ナノ秒単位）。最近のログでのみ利用可能です。\n \n* ArrivalTimeNanosSinceGpsEpoch-GPSエポック（1980/1/6深夜UTC）からのナノ秒の整数。その値は、Rawセンテンスで説明されている一意のエポックごとに、round（（Raw :: TimeNanos --Raw :: FullBiasNanos））に等しくなります。\n \n* RawPseudorangeMeters-メートル単位の生の疑似距離。これは、光速と、信号送信時間（receivedSvTimeInGpsNanos）から信号到着時間（Raw :: TimeNanos-Raw :: FullBiasNanos-Raw ;;biasNanos）までの時間差の積です。その不確実性は、光速とReceivedSvTimeUncertaintyNanosの積で概算できます。\n \n* SignalType-GNSS信号タイプは、コンスタレーション名と周波数帯域の組み合わせです。スマートフォンで測定される一般的な信号タイプには、GPS_L1、GPS_L5、GAL_E1、GAL_E5A、GLO_G1、BDS_B1I、BDS_B1C、BDS_B2A、QZS_J1、QZS_J5があります。\n \n* ReceivedSvTimeNanosSinceGpsEpoch-チップセットが受信した信号伝送時間（GPSエポック以降のナノ秒数）。ReceivedSvTimeNanosから変換されたこの派生値は、すべてのコンステレーションの統一されたタイムスケールにあり、ReceivedSvTimeNanosは、GLONASSの時刻と非GLONASSコンステレーションの時刻を示します。\n \n* SvPosition[X/Y/Z]EcefMeters-ttx = receiveSvTimeInGpsNanosとして定義される「真の信号伝送時間」の最良の推定値でのECEF座標フレーム内の衛星位置（メートル）-satClkBiasNanos（以下に定義）。それらは衛星放送の天体暦で計算され、実際の衛星位置に対して約1メートルの誤差があります。\n \n* Sv[Elevation/Azimuth]Degrees-衛星の高度と方位角（度単位）。これらは、WLSの推定ユーザー位置を使用して計算されます。\n \n* SvVelocity[X/Y/Z]EcefMetersPerSecond-「真の信号伝送時間」ttxの最良の推定値でのECEF座標フレーム内の衛星速度（メートル/秒）。これらは、このアルゴリズムを使用して、衛星放送の天体暦で計算されます。\n \n* SvClockBiasMeters-信号送信時間（receivedSvTimeInGpsNanos）でのメートル単位の衛星ハードウェア遅延と組み合わせた衛星時間補正。その時間に相当するものは、satClkBiasNanosと呼ばれます。satClkBiasNanosは、satelliteTimeCorrectionからsatelliteHardwareDelayを引いたものに等しくなります。IS-GPS-200Hセクション20.3.3.3.3.1で定義されているように、satelliteTimeCorrectionはΔtsv= af0 + af1（t-toc）+ af2（t-toc）2 + Δtrから計算されますが、satelliteHardwareDelayはセクション20.3で定義されています。 3.3.3.2。上記の式のパラメータは、衛星放送の天体暦で提供されています。\n \n* SvClockDriftMetersPerSecond-衛星クロックは、信号送信時間（receivedSvTimeInGpsNanos）でメートル/秒でドリフトします。これは、t+0.5sとt-0.5sでの衛星クロックバイアスの差に等しくなります。\n \n* IsrbMeters-非GPS-L1信号からGPS-L1信号までのメートル単位の信号間範囲バイアス（ISRB）。たとえば、GPS L5のisrbMが1000mの場合、GPSL5疑似距離は同じGPS衛星によって送信されたGPSL1疑似距離よりも1000m長いことを意味します。GPS-L1信号の場合はゼロです。ISRBはGPSチップセットレベルで導入され、加重最小二乗エンジンの状態として推定されます。\n \n* IonosphericDelayMeters-クロブシャーモデルで推定された、メートル単位のイオノスフィア遅延。\n \n* TroposphericDelayMeters-Nigel Penna、Alan Dodson、およびW. Chen（2001）によるEGNOSモデルで推定された、メートル単位の対流圏遅延。\n \n* WlsPositionXEcefMeters--WlsPositionYEcefMeters、WlsPositionZEcefMeters：加重最小二乗（WLS）ソルバーによって推定されたECEF内のユーザー位置。","metadata":{"id":"yexdKx0RKxFy"}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pylab as plt","metadata":{"id":"ekz9bLTyK0Bc","execution":{"iopub.status.busy":"2022-05-09T13:19:36.78924Z","iopub.execute_input":"2022-05-09T13:19:36.789644Z","iopub.status.idle":"2022-05-09T13:19:36.817981Z","shell.execute_reply.started":"2022-05-09T13:19:36.789527Z","shell.execute_reply":"2022-05-09T13:19:36.817236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_device_gnss=pd.read_csv('../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_gnss.csv')\ndf_device_gnss.head()","metadata":{"id":"9O5Mj7gbLCfF","execution":{"iopub.status.busy":"2022-05-09T13:19:36.821025Z","iopub.execute_input":"2022-05-09T13:19:36.822797Z","iopub.status.idle":"2022-05-09T13:19:38.901313Z","shell.execute_reply.started":"2022-05-09T13:19:36.82275Z","shell.execute_reply":"2022-05-09T13:19:38.900335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_device_gnss)","metadata":{"id":"Bn2kg5_lLDtc","execution":{"iopub.status.busy":"2022-05-09T13:19:38.902793Z","iopub.execute_input":"2022-05-09T13:19:38.903096Z","iopub.status.idle":"2022-05-09T13:19:38.909786Z","shell.execute_reply.started":"2022-05-09T13:19:38.903053Z","shell.execute_reply":"2022-05-09T13:19:38.908905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_device_gnss[\"HardwareClockDiscontinuityCount\"].unique()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:42:40.062275Z","iopub.execute_input":"2022-05-09T13:42:40.062565Z","iopub.status.idle":"2022-05-09T13:42:40.070878Z","shell.execute_reply.started":"2022-05-09T13:42:40.062523Z","shell.execute_reply":"2022-05-09T13:42:40.070124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_device_gnss[\"IsrbMeters\"].mean()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:35:10.438432Z","iopub.execute_input":"2022-05-09T13:35:10.439336Z","iopub.status.idle":"2022-05-09T13:35:10.447969Z","shell.execute_reply.started":"2022-05-09T13:35:10.439289Z","shell.execute_reply":"2022-05-09T13:35:10.447137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[train / test] / [drive_id] / [phone_name] /device_imu.csv-電話の加速度計、ジャイロスコープ、磁力計を読み取ります。\n\n* MessageType-行のデータが3つの機器のどれからのものか。\n\n* utcTimeMillis-以下の合計とelapsedRealtimeNanos、最近のNTP（Network Time Protocol）同期後のUTCでの推定デバイス起動時間。\n\n* Measurement[X/Y/Z]-[x / y / z]_uncalib（バイアス補正なし）。\n\n* Bias[X/Y/Z]MicroT-推定[x/y /z]_bias。以前の日付で収集されたデータセットではヌル。","metadata":{}},{"cell_type":"code","source":"df_device_imu=pd.read_csv('../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_imu.csv')\ndf_device_imu.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:38.911872Z","iopub.execute_input":"2022-05-09T13:19:38.912104Z","iopub.status.idle":"2022-05-09T13:19:40.261322Z","shell.execute_reply.started":"2022-05-09T13:19:38.912075Z","shell.execute_reply":"2022-05-09T13:19:40.260466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_device_imu)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:40.262493Z","iopub.execute_input":"2022-05-09T13:19:40.262763Z","iopub.status.idle":"2022-05-09T13:19:40.269713Z","shell.execute_reply.started":"2022-05-09T13:19:40.262733Z","shell.execute_reply":"2022-05-09T13:19:40.268824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_device_imu[\"MessageType\"].unique()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:41:07.80296Z","iopub.execute_input":"2022-05-09T13:41:07.803272Z","iopub.status.idle":"2022-05-09T13:41:07.871494Z","shell.execute_reply.started":"2022-05-09T13:41:07.80324Z","shell.execute_reply":"2022-05-09T13:41:07.870423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(8, 8))\nplt.rcParams[\"font.size\"] = 20\nplt.bar(\n    df_device_imu[\"MessageType\"].value_counts().sort_values(ascending=False).index,\n    df_device_imu[\"MessageType\"].value_counts().sort_values(ascending=False),\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:40:19.203374Z","iopub.execute_input":"2022-05-09T13:40:19.203733Z","iopub.status.idle":"2022-05-09T13:40:19.580393Z","shell.execute_reply.started":"2022-05-09T13:40:19.203694Z","shell.execute_reply":"2022-05-09T13:40:19.579624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"train / [drive_id] /[phone_name]/ground_truth.csv-予想されるタイムスタンプで場所を参照します。\n\n* MessageType-文の接頭辞「Fix」。\n\n* Provider-「GT」、グラウンドトゥルースの略。\n\n* [Latitude/Longitude]Degrees- 参照GNSS受信機（NovAtel SPAN）によって推定されたWGS84の緯度、経度（10進角）。NMEAファイルから抽出する場合、線形補間が適用され、場所が予想される非整数のタイムスタンプに揃えられます。\n\n* AltitudeMeters-参照GNSS受信機によって推定されたWGS84楕円体の上の高さ（メートル単位）。\n\n* SpeedMps*-地上の速度（メートル/秒）。\n\n* AccuracyMeters-68パーセンタイル信頼水準でのこの場所のメートル単位の推定水平精度半径。これは、デバイスの実際の場所が、報告された場所のこの不確実性の距離内にある可能性が68％あることを意味します。\n\n* BearingDegrees-方位は、北から時計回りに度で測定されます。範囲は0〜359.999度です。\n\n* UnixTimeMillis-GPSエポック（1970/1/1真夜中UTC）からのミリ秒の整数。GnssClockから変換されました。","metadata":{}},{"cell_type":"code","source":"df_ground_truth=pd.read_csv('../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/ground_truth.csv')\ndf_ground_truth.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:40.270957Z","iopub.execute_input":"2022-05-09T13:19:40.271196Z","iopub.status.idle":"2022-05-09T13:19:40.307617Z","shell.execute_reply.started":"2022-05-09T13:19:40.271167Z","shell.execute_reply":"2022-05-09T13:19:40.306659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_ground_truth)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:31:11.624993Z","iopub.execute_input":"2022-05-09T13:31:11.625657Z","iopub.status.idle":"2022-05-09T13:31:11.633709Z","shell.execute_reply.started":"2022-05-09T13:31:11.625617Z","shell.execute_reply":"2022-05-09T13:31:11.632784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**霊夢:地図上に通った道筋を書いてみよう。Rob氏のコードを参考にするよ。**  \n\n**Reimu: Let's write the path you took on the map. I'll refer to Rob's code.**\n\nhttps://www.kaggle.com/code/robikscube/smartphone-competition-2022-twitch-stream","metadata":{}},{"cell_type":"code","source":"!pip install nb_black > /dev/null","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:40.309647Z","iopub.execute_input":"2022-05-09T13:19:40.309909Z","iopub.status.idle":"2022-05-09T13:19:55.536255Z","shell.execute_reply.started":"2022-05-09T13:19:40.309868Z","shell.execute_reply":"2022-05-09T13:19:55.535015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%load_ext lab_black","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:55.54032Z","iopub.execute_input":"2022-05-09T13:19:55.540622Z","iopub.status.idle":"2022-05-09T13:19:55.860859Z","shell.execute_reply.started":"2022-05-09T13:19:55.540586Z","shell.execute_reply":"2022-05-09T13:19:55.859751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nimport glob\nfrom dataclasses import dataclass\nfrom tqdm.notebook import tqdm\nfrom scipy.interpolate import InterpolatedUnivariateSpline\n\npd.set_option(\"max_columns\", 500)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:55.863472Z","iopub.execute_input":"2022-05-09T13:19:55.863739Z","iopub.status.idle":"2022-05-09T13:19:57.611275Z","shell.execute_reply.started":"2022-05-09T13:19:55.863708Z","shell.execute_reply":"2022-05-09T13:19:57.610312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/ground_truth.csv\"\n)\ngnss = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_gnss.csv\"\n)\nimu = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_imu.csv\"\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:57.61511Z","iopub.execute_input":"2022-05-09T13:19:57.615348Z","iopub.status.idle":"2022-05-09T13:19:59.651858Z","shell.execute_reply.started":"2022-05-09T13:19:57.615319Z","shell.execute_reply":"2022-05-09T13:19:59.650775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_PATH = \"../input/smartphone-decimeter-2022\"\n\nWGS84_SEMI_MAJOR_AXIS = 6378137.0\nWGS84_SEMI_MINOR_AXIS = 6356752.314245\nWGS84_SQUARED_FIRST_ECCENTRICITY = 6.69437999013e-3\nWGS84_SQUARED_SECOND_ECCENTRICITY = 6.73949674226e-3\n\nHAVERSINE_RADIUS = 6_371_000\n\n\n@dataclass\nclass ECEF:\n    x: np.array\n    y: np.array\n    z: np.array\n\n    def to_numpy(self):\n        return np.stack([self.x, self.y, self.z], axis=0)\n\n    @staticmethod\n    def from_numpy(pos):\n        x, y, z = [np.squeeze(w) for w in np.split(pos, 3, axis=-1)]\n        return ECEF(x=x, y=y, z=z)\n\n\n@dataclass\nclass BLH:\n    lat: np.array\n    lng: np.array\n    hgt: np.array\n\n\ndef ECEF_to_BLH(ecef):\n    a = WGS84_SEMI_MAJOR_AXIS\n    b = WGS84_SEMI_MINOR_AXIS\n    e2 = WGS84_SQUARED_FIRST_ECCENTRICITY\n    e2_ = WGS84_SQUARED_SECOND_ECCENTRICITY\n    x = ecef.x\n    y = ecef.y\n    z = ecef.z\n    r = np.sqrt(x**2 + y**2)\n    t = np.arctan2(z * (a / b), r)\n    B = np.arctan2(z + (e2_ * b) * np.sin(t) ** 3, r - (e2 * a) * np.cos(t) ** 3)\n    L = np.arctan2(y, x)\n    n = a / np.sqrt(1 - e2 * np.sin(B) ** 2)\n    H = (r / np.cos(B)) - n\n    return BLH(lat=B, lng=L, hgt=H)\n\n\ndef haversine_distance(blh_1, blh_2):\n    dlat = blh_2.lat - blh_1.lat\n    dlng = blh_2.lng - blh_1.lng\n    a = (\n        np.sin(dlat / 2) ** 2\n        + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng / 2) ** 2\n    )\n    dist = 2 * HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\n\ndef pandas_haversine_distance(df1, df2):\n    blh1 = BLH(\n        lat=np.deg2rad(df1[\"LatitudeDegrees\"].to_numpy()),\n        lng=np.deg2rad(df1[\"LongitudeDegrees\"].to_numpy()),\n        hgt=0,\n    )\n    blh2 = BLH(\n        lat=np.deg2rad(df2[\"LatitudeDegrees\"].to_numpy()),\n        lng=np.deg2rad(df2[\"LongitudeDegrees\"].to_numpy()),\n        hgt=0,\n    )\n    return haversine_distance(blh1, blh2)\n\n\ndef ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis):\n    ecef_columns = [\n        \"WlsPositionXEcefMeters\",\n        \"WlsPositionYEcefMeters\",\n        \"WlsPositionZEcefMeters\",\n    ]\n    columns = [\"utcTimeMillis\"] + ecef_columns\n    ecef_df = (\n        gnss_df.drop_duplicates(subset=\"utcTimeMillis\")[columns]\n        .dropna()\n        .reset_index(drop=True)\n    )\n    ecef = ECEF.from_numpy(ecef_df[ecef_columns].to_numpy())\n    blh = ECEF_to_BLH(ecef)\n\n    TIME = ecef_df[\"utcTimeMillis\"].to_numpy()\n    lat = InterpolatedUnivariateSpline(TIME, blh.lat, ext=3)(UnixTimeMillis)\n    lng = InterpolatedUnivariateSpline(TIME, blh.lng, ext=3)(UnixTimeMillis)\n    return pd.DataFrame(\n        {\n            \"tripId\": tripID,\n            \"UnixTimeMillis\": UnixTimeMillis,\n            \"LatitudeDegrees\": np.degrees(lat),\n            \"LongitudeDegrees\": np.degrees(lng),\n        }\n    )\n\n\ndef calc_score(tripID, pred_df, gt_df):\n    d = pandas_haversine_distance(pred_df, gt_df)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])\n    return score","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:59.653346Z","iopub.execute_input":"2022-05-09T13:19:59.653614Z","iopub.status.idle":"2022-05-09T13:19:59.755882Z","shell.execute_reply.started":"2022-05-09T13:19:59.653582Z","shell.execute_reply":"2022-05-09T13:19:59.754848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_traffic(\n    df,\n    lat_col=\"LatitudeDegrees\",\n    lon_col=\"LongitudeDegrees\",\n    center=None,\n    color_col=\"phone\",\n    label_col=\"tripId\",\n    zoom=9,\n    opacity=1,\n):\n    if center is None:\n        center = {\n            \"lat\": df[lat_col].mean(),\n            \"lon\": df[lon_col].mean(),\n        }\n    fig = px.scatter_mapbox(\n        df,\n        lat=lat_col,\n        lon=lon_col,\n        color=color_col,\n        labels=label_col,\n        zoom=zoom,\n        center=center,\n        height=600,\n        width=800,\n        opacity=0.5,\n    )\n    fig.update_layout(mapbox_style=\"stamen-terrain\")\n    fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n    fig.update_layout(title_text=\"GPS trafic\")\n    fig.show()\n\n\ndef plot_gt_vs_baseline(tripId):\n    \"\"\"\n    Create a plot of the baseline predictions vs. the ground truth\n    for a given tripId\n    \"\"\"\n    # Pull Data for an example phone\n    gt = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/ground_truth.csv\"\n    )\n    gnss = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/device_gnss.csv\"\n    )\n    imu = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/device_imu.csv\"\n    )\n    baseline = ecef_to_lat_lng(trip_id, gnss, gt[\"UnixTimeMillis\"].values)\n    # Combine ground truth with baseline predictions\n    baseline[\"isGT\"] = False\n    gt[\"isGT\"] = True\n    gt[\"tripId\"] = tripId\n\n    combined = (\n        pd.concat([baseline, gt[baseline.columns]], axis=0)\n        .reset_index(drop=True)\n        .copy()\n    )\n\n    # Plotting the route\n    visualize_traffic(\n        combined,\n        lat_col=\"LatitudeDegrees\",\n        lon_col=\"LongitudeDegrees\",\n        color_col=\"isGT\",\n        zoom=10,\n    )","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:59.757174Z","iopub.execute_input":"2022-05-09T13:19:59.75742Z","iopub.status.idle":"2022-05-09T13:19:59.806708Z","shell.execute_reply.started":"2022-05-09T13:19:59.75739Z","shell.execute_reply":"2022-05-09T13:19:59.805803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trip_id = \"2020-05-15-US-MTV-1/GooglePixel4XL\"\nplot_gt_vs_baseline(trip_id)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:19:59.807893Z","iopub.execute_input":"2022-05-09T13:19:59.808102Z","iopub.status.idle":"2022-05-09T13:20:03.291982Z","shell.execute_reply.started":"2022-05-09T13:19:59.808076Z","shell.execute_reply":"2022-05-09T13:20:03.291082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trip_id2 = \"2020-05-21-US-MTV-1/GooglePixel4\"\nplot_gt_vs_baseline(trip_id2)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:03.293366Z","iopub.execute_input":"2022-05-09T13:20:03.293718Z","iopub.status.idle":"2022-05-09T13:20:05.519378Z","shell.execute_reply.started":"2022-05-09T13:20:03.293681Z","shell.execute_reply":"2022-05-09T13:20:05.51827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trip_id3 = \"2020-08-06-US-MTV-1/GooglePixel4\"\nplot_gt_vs_baseline(trip_id3)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:05.520639Z","iopub.execute_input":"2022-05-09T13:20:05.520899Z","iopub.status.idle":"2022-05-09T13:20:07.168783Z","shell.execute_reply.started":"2022-05-09T13:20:05.520868Z","shell.execute_reply":"2022-05-09T13:20:07.167297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# metadata","metadata":{}},{"cell_type":"code","source":"df_meta = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/metadata/constellation_type_mapping.csv\"\n)\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:07.170384Z","iopub.execute_input":"2022-05-09T13:20:07.170697Z","iopub.status.idle":"2022-05-09T13:20:07.201632Z","shell.execute_reply.started":"2022-05-09T13:20:07.170663Z","shell.execute_reply":"2022-05-09T13:20:07.200428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# test","metadata":{}},{"cell_type":"markdown","source":"**霊夢:次はTestを見てみよう**  \n\n**Reimu: Let's take a look at Test next**","metadata":{"id":"kQzI6s13NdgP"}},{"cell_type":"code","source":"df_test_gnss = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/test/2021-04-28-US-MTV-2/SamsungGalaxyS20Ultra/device_gnss.csv\"\n)\ndf_test_gnss.head()","metadata":{"id":"KyHHd-5rMFH9","execution":{"iopub.status.busy":"2022-05-09T13:20:07.203228Z","iopub.execute_input":"2022-05-09T13:20:07.203476Z","iopub.status.idle":"2022-05-09T13:20:08.423719Z","shell.execute_reply.started":"2022-05-09T13:20:07.203446Z","shell.execute_reply":"2022-05-09T13:20:08.422478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_test_gnss)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:08.425412Z","iopub.execute_input":"2022-05-09T13:20:08.425779Z","iopub.status.idle":"2022-05-09T13:20:08.436129Z","shell.execute_reply.started":"2022-05-09T13:20:08.425734Z","shell.execute_reply":"2022-05-09T13:20:08.435045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test_gnss[\"SignalType\"].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:29:48.902365Z","iopub.execute_input":"2022-05-09T13:29:48.902699Z","iopub.status.idle":"2022-05-09T13:29:48.921955Z","shell.execute_reply.started":"2022-05-09T13:29:48.902649Z","shell.execute_reply":"2022-05-09T13:29:48.921088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test_gnss[\"CodeType\"].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:50:16.376418Z","iopub.execute_input":"2022-05-09T13:50:16.37674Z","iopub.status.idle":"2022-05-09T13:50:16.394514Z","shell.execute_reply.started":"2022-05-09T13:50:16.376709Z","shell.execute_reply":"2022-05-09T13:50:16.393823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test_imu = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/test/2021-04-28-US-MTV-2/SamsungGalaxyS20Ultra/device_imu.csv\"\n)\ndf_test_imu.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:08.437783Z","iopub.execute_input":"2022-05-09T13:20:08.438117Z","iopub.status.idle":"2022-05-09T13:20:09.302606Z","shell.execute_reply.started":"2022-05-09T13:20:08.438073Z","shell.execute_reply":"2022-05-09T13:20:09.300908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(df_test_imu)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T13:20:09.304875Z","iopub.execute_input":"2022-05-09T13:20:09.305513Z","iopub.status.idle":"2022-05-09T13:20:09.313966Z","shell.execute_reply.started":"2022-05-09T13:20:09.30546Z","shell.execute_reply":"2022-05-09T13:20:09.313094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# sample_submission","metadata":{}},{"cell_type":"markdown","source":"**魔理沙:最後にsample_submissionを見てみよう**  \n\n**Marisa: Finally, let's take a look at sumple_submission**","metadata":{"id":"whxCx7WAN_us"}},{"cell_type":"code","source":"sumple_submission = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/sample_submission.csv\"\n)\nsumple_submission.head()","metadata":{"id":"3f2UQrtJNYFg","execution":{"iopub.status.busy":"2022-05-08T13:41:54.445497Z","iopub.execute_input":"2022-05-08T13:41:54.445997Z","iopub.status.idle":"2022-05-08T13:41:54.519193Z","shell.execute_reply.started":"2022-05-08T13:41:54.445964Z","shell.execute_reply":"2022-05-08T13:41:54.518376Z"},"trusted":true},"execution_count":null,"outputs":[]}]}