7 地理データの再投影

必須パッケージ

  • この章では、以下のパッケージが必要である。

7.1 イントロダクション

Section 2.4 は、座標参照系 (CRS) を紹介し、2つの主要なタイプに焦点を当てた。地理座標系(‘lon/lat’、単位は経度と緯度)と投影座標系(通常は基準点からのメートル)である。 この章では、その知識をもとに、さらに踏み込んだ内容になっている。 本書では、ある CRS から別の CRS に地理データを設定し、変換する方法を説明する。さらに、特にデータが緯度経度座標で保存されている場合、CRS を無視することによって発生し得る問題があるので、これを明らかにして注意を促したい。

多くのプロジェクトでは、CRS について心配する必要はないし、変換について考える必要もない。 しかし、データが投影座標系なのか地理座標系なのか、そしてそれがジオメトリ操作に与える影響を知ることは重要である。 データの CRS とジオメトリ操作の結果(次のセクションで取り上げる)を知っていれば、CRS は裏でただうまく機能してくれる。しかし、物事がうまくいかないとき、突然 CRS について学ぶ必要が出てくるのである。 すべての投影データが入る投影 CRS を明確に定義し、異なる CRS をどのように、そしてなぜ使うのかを理解することで、物事がうまくいくことを確実にすることができる。 さらに、座標系について学ぶことで、地理データセットとその効果的な使用方法についての知識を深めることができる。

この章では、CRS の基本を学び、異なる CRS を使用した場合の結果(何が問題になるかを含む)を示し、ある座標系から別の座標系にデータセットを「再投影」する方法について説明する。 次のセクションでは R における CRS を紹介し、続いて Section 7.2 で空間オブジェクトに関連する CRS の取得と設定方法を示す。 Section 7.4 は、バッファを作成する作業例を参照しながら、データがどの CRS にあるのかを知ることの重要性を示している。 Section 7.5 と Section 7.6 において、それぞれ、いつ再投影するべきか、どの CRS を使うかという問題に取り組んでいる。 ベクタとラスタの再投影については、7.77.8 で、地図投影の修正については Section 7.9 で説明している。

7.2 座標参照系

R-spatial のコアパッケージや QGIS などのデスクトップ GIS ソフトウェアなど、CRS 変換を必要とする最新の地理ツールのほとんどは、「ある座標参照系(CRS)から別の座標に変換する」オープンソース C++ ライブラリ PROJ とつながっている。 CRS には、以下のような様々な表現がある。

  1. 単純だが曖昧になる可能性のある記述、例えば「lon/lat座標で表示される」。
  2. +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs のような、形式化された、しかし今では時代遅れの「proj4 文字列」。
  3. EPSG:4326 のような識別用の ‘authority:code’ テキスト文字列。

上の例はそれぞれ、全地球測位システム(GPS)座標やその他多くのデータセットの基礎となる「WGS84」座標系という同じものを指している。(訳注:WGS84 は、世界で一般的に用いられている緯度経度座標系だが、日本の場合は、JGD2011 が用いられることが多い。そのコードは EPSG:6668 である。なお、JGD2011 は東日本大震災の影響を考慮しており、震災以前は JGD2000 が用いられていた。そのコードは EPSG:4612 である。データが作成された年によって両者を使い分けるとよい。詳細については、空間情報クラブなどを参照。) しかし、どれが正しいのだろうか?

短く答えると、CRS を識別する3つ目の方法が正しい。EPSG:4326 は、本書で取り上げた sfterra (さらに stars)パッケージ、及び QGISPROJ などの地理データを扱うための多くのソフトウェアプロジェクトによって理解されている。 EPSG:4326 は将来を見据えたものである。 さらに、機械可読でありながら、proj-string 表現とは異なり、「EPSG:4326」は短く、覚えやすく、オンラインで非常に「見つけやすい」(例えば、EPSG:4326 を検索すると、ウェブサイト epsg.io の専用ページが表示される)。 より簡潔に 4326 だけでも sf によって理解されるが、 曖昧さを防ぎ、文脈を提供するために、より明示的な AUTHORITY:CODE という形式を推奨する

より長く答えると、3つの記述のどれも十分ではなく、CRS の処理と変換を明確にするためには、より詳細な情報が必要になる。 このため、Open Geospatial Consortium(OGC、sf パッケージが実装するシンプルフィーチャの仕様も整備している団体)が、WKT(Well Known Text)と呼ばれる CRS 記述形式をオープンスタンダードで開発した。 これは、「ISO 19111:2019 に記述された座標参照系の抽象モデルのテキスト文字列実装の構造と内容を定義する」100ページを超えるドキュメントに詳細が記述されている。 (Open Geospatial Consortium 2019)。 WGS84 CRS、識別子 EPSG:4326 の WKT 表現は以下の通りである。

st_crs("EPSG:4326")
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["Horizontal component of 3D system."],
#>         AREA["World."],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]

コマンドの出力は、CRS 識別子(空間参照識別子または SRIDとも呼ばれる)がどのように機能するかを示している。これは単にルックアップであり、CRS のより完全な WKT 表現に関連する一意の識別子を提供するものである。 このことは、識別子と CRS の長い WKT 表現との間にミスマッチがある場合はどうなるのか、という問題を提起している。 この点、Open Geospatial Consortium (2019) は明確で、冗長な WKT 表現が識別子より優先される。

引用された識別子の属性や値が、WKT 記述で明示的に与えられた属性や値と矛盾する場合、 WKT の値が優先されるものとする。

CRS の識別子を AUTHORITY:CODE という形式で参照する慣習は、他の言語で書かれた地理ソフトウェアでも使われており、正式に定義された広範囲の座標系を参照することができる。26 CRS の識別子で最もよく使われる機関は、CRS の標準化リストを発表した欧州石油調査グループ(European Petroleum Survey Group)の頭文字をとった EPSG である(EPSG は2005年に石油・ガス団体の Geomatics Committee of the International Association of Oil & Gas Producers により引き継ぎされた)。 CRS の識別子には、他の機関を使用することもできる。 例えば、ESRI:54030 は、ESRI の Robinson 投影の実装で、以下の WKT 文字列(最初の8行のみ表示)を持っていることを指している。

sf::st_crs("ESRI:54030")
#> Coordinate Reference System:
#>   User input: ESRI:54030 
#>   wkt:
#> PROJCRS["World_Robinson",
#>     BASEGEOGCRS["WGS 84",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]]],
#> ...

WKT 文字列は網羅的で詳細かつ正確であるため、CRS の格納や変換を曖昧にすることなく行うことができる。 測地系、楕円体、本初子午線、投影法、単位など、任意の CRS に関するすべての関連情報が含まれている。27

最近の PROJ のバージョン(6+)では、座標操作を定義するために proj-string を使用することができるが、いくつかの proj-string キー(+nadgrids , +towgs84 , +k , +init=epsg:)はもうサポートされていないか、推奨されないものである。 また、proj-string に直接設定できる測地系は3つ(WGS84、NAD83、NAD27)だけである。 CRS の定義の進化と PROJ ライブラリの長い説明は、Bivand (2021)Pebesma and Bivand (2022) の第2章、Floris Vanderhaeghe によるブログに記載されている。 PROJ documentation に概説されているように、WKT CRS 形式には WKT1 と 2 種類の WKT2 があり、後者 (WKT2, 2018 仕様) は ISO 19111:2019 に対応するものである (Open Geospatial Consortium 2019)

7.3 座標系の照会と設定

R の空間オブジェクトに CRS がどのように格納され、どのようにクエリや設定ができるのかを見ていこう。 まず、ベクタの地理データオブジェクトの CRS の取得と設定について、次の例から見ていく。

vector_filepath = system.file("shapes/world.gpkg", package = "spData")
new_vector = read_sf(vector_filepath)

新しいオブジェクトである new_vector は、世界の国々を表すクラス sf のデータフレームである(詳しくはヘルプページ ?spData::world を参照)。 CRS は sf 関数 st_crs() で取得することができる。

st_crs(new_vector) # get CRS
#> Coordinate Reference System:
#>   User input: WGS 84 
#>   wkt:
#>   ...

出力は主に2つの要素を含むリストである。

  1. User input (この場合、WGS 84、入力ファイルから取得した EPSG:4326 の同義語)、前述の CRS 識別子に対応する。
  2. wkt CRS に関するすべての関連情報を含む完全な WKT 文字列を含む。

input 要素は柔軟性があり、入力ファイルやユーザー入力に応じて、AUTHORITY:CODE 表現(例:EPSG:4326)、CRS の名前(例: WGS 84)、あるいは proj-string 定義を含めることができる。 wkt 要素には WKT 表現が格納され、オブジェクトをファイルに保存したり、座標演算を行う際に使用される。 上記で、new_vector のオブジェクトは、WGS84 楕円体を持ち、グリニッジ本初子午線を使用し、緯度・経度軸の順番になっていることがわかる。 この場合、この CRS の使用に適したエリアを説明する USAGE や、CRS の識別子を指し示す IDEPSG:4326) などの追加要素もある。

また、st_crs 関数には、使用中の CRS に関する追加情報を取得することができる、という便利な特徴もある。 例えば、実行させてみよう。

  • st_crs(new_vector)$IsGeographic、CRS がジオグラフィックかどうかを確認する。
  • st_crs(new_vector)$units_gdal、CRS 単位を調べることができる。
  • st_crs(new_vector)$srid その「SRID」識別子を抽出する(利用可能な場合)。
  • st_crs(new_vector)$proj4string proj-string表現を抽出する。

座標参照系(CRS)がない場合や間違った CRS が設定されている場合は、st_set_crs() 関数を使用することができる(この場合、ファイル読み込み時にすでに CRS が正しく設定されているので、WKT 文字列は変更されずに残る)。

new_vector = st_set_crs(new_vector, "EPSG:4326") # CRS を設定

CRS の取得と設定は、ラスタ 地理データオブジェクトと同様の方法で行われる。 terra パッケージの crs() 関数は SpatRaster オブジェクトから CRS 情報にアクセスする (きれいに印刷するために cat() 関数を使用していることに注意)。

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
cat(crs(my_rast)) # CRS を取得
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     ID["EPSG",4326]]

出力は CRS の WKT 表現である。 同じ関数、crs() を使って、ラスタオブジェクトに CRS を設定することもできる。

crs(my_rast) = "EPSG:26912" # CRS を設定

ここでは、識別子(ほとんどの場合推奨)または完全な WKT 表現のいずれかを使用することができる。 crs を設定する代替方法としては、proj-string 文字列または crs() を持つ他の既存のオブジェクトから抽出された CRS があるが、これらのアプローチは将来対応されない可能性がある。

重要なのは、st_crs()crs() 関数は、座標の値や形状を変更しないことである。 その役割は、オブジェクト CRS のメタデータ情報を設定することのみである。

Section 2.2 で紹介したロンドンの例を基に、以下のコードで作成した london データセットのように、地理的オブジェクトの CRS が不明な場合がある。

london = data.frame(lon = -0.1, lat = 51.5) |> 
  st_as_sf(coords = c("lon", "lat"))
st_is_longlat(london)
#> [1] NA

この出力 NA は、sf が CRS が何であるかを知らず、推測するつもりがないことを示している(NA は、“Not Applicable/Available” の略語で、文字通り「利用できない」という意味)。 CRS を手動で指定するか、CRS メタデータを持つソースから読み込まれない限り、sf は座標系について「わからない」と言う以外の明示的な仮定をしない。 この動作は、利用可能な CRS の多様性を考えると理にかなっているが、GeoJSON ファイルフォーマット仕様のような、すべての座標が lon/lat CRS(EPSG:4326`)を持つという単純化した仮定をするいくつかのアプローチとは異なる。

sf オブジェクトに CRS を追加するには、主に 3 つの方法がある。

  • CRS をあらかじめ存在するオブジェクトに割り当てる。例えば st_crs(london) = "EPSG:4326"
  • st_as_sf(... crs = "EPSG:4326") などのジオメトリオブジェクトを作成する sf 関数の crs 引数に CRS を渡すことで、ジオメトリオブジェクトを作成することができる。また、ラスタデータセットを作成する際にも、同じ引数で CRS を設定することができる(例:rast(crs = "EPSG:4326"))。
  • 新しい CRS を持つデータのバージョンを返す st_set_crs()、例を次のコードチャンクで示す。
london_geo = st_set_crs(london, "EPSG:4326")
st_is_longlat(london_geo)
#> [1] TRUE

CRS の指定されていないデータセットは、問題を起こすことがある。すべての地理座標には座標系があり、ソフトウェアでは、使用する CRS の種類がわかっている場合にのみ、プロットやジオメトリ操作に関する適切な判断を行うことができる。

7.4 投影データおよび非投影データに対する幾何学操作

sf version 1.0.0 より、R は緯度経度 CRS を持って入るベクタデータセットに対する機能が大幅に強化された。この機能強化は、Section 2.2.9 で取り上げた S2 球面ジオメトリエンジンによるものである。 Figure 7.1 で示すように、sf は、CRS 種別に応じて GEOS または S2 を使い分ける(デフォルトは S2)。 投影されたデータと CRS がないデータの場合、常に GEOS が使われる。地理データではデフォルトで S2 が使われるが、無効化したい場合は sf::sf_use_s2(FALSE) とする。

入力データの CRS に依存する sf パッケージのジオメトリ操作の動作。

FIGURE 7.1: 入力データの CRS に依存する sf パッケージのジオメトリ操作の動作。

CRS の重要性を示すために、このセクションでは、前セクションで作成した london オブジェクトの周りに 100 km のバッファを作成する。 また、100 km にほぼ相当する1度(赤道では1度は約 111 km)の「距離」を持つ意図的に欠陥のあるバッファを作成する。 コードに入る前に、Figure 7.2 を見て、これからコードチャンクで再現するはずの出力を視覚的に把握するのもよいだろう。

最初の段階として、上記で作成した londonlondon_geo のオブジェクトの周りに、ロンドン中心部から1度と 100 km(または 10万 m、科学的表記法では 1e5 と表現できる)の境界距離の3つのバッファを作成する。

london_buff_no_crs = st_buffer(london, dist = 1)  # incorrect: no CRS
london_buff_s2 = st_buffer(london_geo, dist = 1e5) # silent use of s2
london_buff_s2_100_cells = st_buffer(london_geo, dist = 1e5, max_cells = 100) 

上の1行目では、sf は入力が投影されていると仮定して、度数単位のバッファを持つ結果を生成しているが、後述するようにこれは問題である。 2行目の sf では、Chapter 2 で導入された球面幾何エンジン S2 を無言で使用して、max_cells = 1000 のデフォルト値 — 3行目で 100 に設定 — を使用してバッファの範囲を計算しているが、その結果はすぐに明らかになるだろう(詳しくは ?s2::s2_buffer_cells を参照)。 非投影(地理)座標系に対する sf の S2 幾何エンジン使用の影響を強調するために、以下のコードチャンクで、コマンド sf_use_s2() (デフォルトではオン、TRUE )で一時的にそれを無効にしてみよう。 london_buff_no_crs と同様、新しい london_geo オブジェクトは、単位が度でほとんどの場合意味を持たない。

sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
london_buff_lonlat = st_buffer(london_geo, dist = 1) # incorrect result
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on

上記の警告メッセージは、lon/lat データに対して平面ジオメトリ演算を実行する際の問題を示唆している。 球形幾何演算をコマンド(sf::sf_use_s2(FALSE))でオフにすると、バッファ(およびその他のジオメトリ演算)は緯度と経度の単位を使用するため、メートルなどの適切な距離単位の代わりにならないため、価値のない出力となる場合がある。

子午線と呼ばれる 2 本の経線間の距離は、赤道では約 111 km である(正確な距離は geosphere::distGeo(c(0, 0), c(1, 0)) を実行して調べられる)。 これは極点でゼロに縮まる。 例えば、ロンドンの緯度では、子午線の距離は 70 km 以下である(課題:これを検証するコードを実行しなさい)。 一方、緯度線は緯度に関係なく等距離にあり、赤道や極付近を含めて常に 111 km 程度の間隔がある(Figure 7.2 から Figure 7.4参照)。

地理的(緯度経度)CRS に関する警告は「CRS を設定してはいけない」という解釈をするべきではない。むしろ、ほとんど常に設定すべきである。 これは、投影された CRS にデータを再投影することを提案したと理解するのがよいだろう。 ただし、この提案に常に耳を傾ける必要もない。空間演算と幾何演算を実行しても、ほとんど違いがない場合もある(例えば、空間的な部分集合作成など)。 しかし、バッファ作成など距離を伴う操作では、(球面幾何エンジンを使わずに)良い結果を得るには、データの投影コピーを作成し、それに対して操作を実行するしかない。 これは、以下のコードチャンクで行われる。

london_proj = data.frame(x = 530000, y = 180000) |> 
  st_as_sf(coords = 1:2, crs = "EPSG:27700")

結果は、london と同じであるが、適切な CRS(この場合、EPSG コードが 27700 の British National Grid)に再投影され、単位がメートルである新しいオブジェクトになる。 CRS が変化したことは、st_crs() を使って次のように確認できる(出力の一部は ... で置き換えている)。

st_crs(london_proj)
#> Coordinate Reference System:
#>   User input: EPSG:27700 
#>   wkt:
#> PROJCRS["OSGB36 / British National Grid",
#>     BASEGEOGCRS["OSGB36",
#>         DATUM["Ordnance Survey of Great Britain 1936",
#>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
#>                 LENGTHUNIT["metre",1]]],
#> ...

この CRS の記述で注目すべきは、EPSG コード (EPSG: 27700) と詳細な wkt の文字列(最初の5行のみ表示)であろう。28 LENGTHUNIT フィールドに記述された CRS の単位が(度ではなく)メートルであることから、これが投影型 CRS であることがわかる。 st_is_longlat(london_proj) は現在 FALSE を返し、london_proj に対する幾何演算は警告なしで機能する。 london_proj のバッファ操作は GEOS を使用し、結果は適切な距離の単位で返される。 次のコードは、ちょうど 100 km の投影データの周りにバッファを作成するものである。

london_buff_projected = st_buffer(london_proj, 1e5)

先に作成した CRS を持つ3つの london_buff* オブジェクト(london_buff_s2london_buff_lonlatlondon_buff_projected)の形状を Figure 7.2 に示す。

ロンドン周辺のバッファで、S2球面幾何エンジンを用いて作成した緯度経度データ(左)、投影データ(中)、球面幾何を用いない緯度経度データ(右)の結果を示している。左のプロットは、投影されていないデータを sf でバッファ作成した結果を示しており、デフォルトで Google の S2 spherical geometry engine を max_cells を 1000 に設定して呼び出している(細線)。太い「ブロック状の」線は、max_cells を 100 に設定して同じ操作を行った結果を示している。

FIGURE 7.2: ロンドン周辺のバッファで、S2球面幾何エンジンを用いて作成した緯度経度データ(左)、投影データ(中)、球面幾何を用いない緯度経度データ(右)の結果を示している。左のプロットは、投影されていないデータを sf でバッファ作成した結果を示しており、デフォルトで Google の S2 spherical geometry engine を max_cells を 1000 に設定して呼び出している(細線)。太い「ブロック状の」線は、max_cells を 100 に設定して同じ操作を行った結果を示している。

s2 と適切に予測された CRS に基づくバッファは、バッファの境界のすべての部分がロンドンから等距離にあることを意味し、「つぶされた」ものではないことは、Figure 7.2 から明らかである。 入力に CRS がないか、sf_use_s2() がオフになっているため、s2使われていないときに緯度経度 CRS から生成される結果は、南北軸に細長く歪んでいる。緯度経度に投影データを仮定する(GEOS による)アルゴリズムを使うことが危険であることが明らかである。 しかし、S2 で生成された結果も、劇的な変化はないものの、歪んでいる。 Figure 7.2 (左) のバッファ境界はどちらもギザギザしているが、これは s2 の引数 max_cells を 100 に設定して作成したバッファを表す太い境界の場合のみ、明らかまたは関連性があると考えられる。 S2 経由の緯度経度データから得られる結果は、投影データから得られる結果とは異なるという教訓である。 投影データにおける S2 由来のバッファと GEOS 由来のバッファの差は、max_cells の値が大きくなるほど小さくなる。この引数の「正しい」値は多くの要因に依存すると思われるが、デフォルト値1000は妥当なデフォルト値だと思われる。 max_cells の値を選択する時、多くの場合、計算の速度と結果の解像度のバランスをとる必要がある。 曲線の境界が有利な場合、バッファ作成(または他のジオメトリ操作)の前に投影型 CRS に変換することが適切な場合がある。

CRS の重要性(主に投影か地理的か)と、緯度経度に対してバッファ作成する際に S2 を使用するという sf のデフォルト設定がデータに与える影響は、上記の例から明らかである。 この後のセクションでは、投影 CRS が必要な場合にどの CRS を使用するか、ベクタおよびラスタオブジェクトの再投影の詳細について、より深く掘り下げて説明する。

7.5 いつ再投影するべきか?

前節では、CRS を手動で設定する方法として、st_set_crs(london, "EPSG:4326") を紹介した。 しかし、現実のアプリケーションでは、データの読み込み時に自動的にCRS が設定されるのが一般的である。 多くのプロジェクトで、CRS 関連の主なタスクは、ある CRS から別の CRS に、オブジェクトを変換することである。 しかし、どのような場合にデータを変換する必要があるのだろうか。 そして、どの CRS に? こういった質問に対する明確な答えはなく、CRS の選択には常にメリットだけでなくデメリットもある (Maling 1992)。 このセクションでは、決定する際に役立つ一般原則を紹介しよう。

まず最初に、いつ変換するべきかを考える。 Figure 7.2 が示すように、st_buffer() など幾何学的な関数を使う場合、投影型 CRS への変換が不可欠な場合がある。 逆に、leaflet パッケージでデータをオンライン公開する場合は、地理的 CRS が必要になる。 また、異なる CRS を持つ2つのオブジェクトの距離を求める場合のように、異なる CRS を持つ2つのオブジェクトを比較したり、組み合わせたりする必要がある場合もある。

st_distance(london_geo, london_proj)
# > Error: st_crs(x) == st_crs(y) is not TRUE

londonlondon_proj のオブジェクトを地理的に比較できるようにするには、一方を他方の CRS に変換する必要がある。 しかし、どの CRS を使えばいいのか? 特にウェブマッピングを含む多くのプロジェクトでは、EPSG:4326 での出力が必要であるが、その場合、投影オブジェクトを変換する価値がある。 しかし、球面幾何演算エンジンではなく平面ジオメトリ演算が必要なプロジェクト(例えば、滑らかなエッジを持つバッファを作成する)の場合、地理的 CRS のデータを英国ナショナルグリッド(EPSG:27700)などの投影 CRS で同等のオブジェクトに変換する価値があるだろう。 それが、Section 7.7 のテーマである。

7.6 どの CRS を使うべきか?

どのCRSを使うかというのは難しい問題で、「正しい」答えがあるわけではない。 「万能の投影は存在せず、すべて指定したフレームの中心から離れると歪みが発生する」 (Bivand, Pebesma, and Gómez-Rubio 2013) 。 さらに言えば、すべてのタスクで1つのプロジェクターだけに執着するべきではない。 ある投影法を解析の一部に使い、別の投影法を別の部分に使い、さらに別の投影法を可視化することも可能である。 常に自分の目標に最も適した CRS を選ぶように心がけよう

地理的 CRS を選択する場合、WGS84 という答えになることが多いようである。(訳注:上で既に述べたが、WGS84 に対応する日本の CRS は、東日本大震災以降の JGD2011 と、それ以前の JGD2000 である。) ウェブマッピングだけでなく、GPS データセットや何千ものラスタ、ベクタデータセットがこの CRS でデフォルトで提供されているため、利用されている。 WGS84 は世界で最も一般的な CRS なので、その EPSG コード 4326 を知っておくとよいだろう。(訳注:JGD2011は 6668、JGD2000 は 4612 である。) この「マジックナンバー」は、投影された CRS が異常なオブジェクトを、広く理解されるものに変換するために使用することができる。

投影 CRS が必要な場合はどうだろうか? 場合によっては、自由に決められるものではない。 「多くの場合、投影の選択は公的な地図作成機関によって行われる」 (Bivand, Pebesma, and Gómez-Rubio 2013) 。 つまり、現地のデータソースで作業する場合、公式の CRS が最も正確ではないとしても、互換性を確保するために、データが提供された CRS で作業することが望ましいと思われる。 ロンドンの例は、(a) 英国ナショナルグリッド(関連する EPSG コード27700)がよく知られており、(b) 元のデータセット(london)がすでにその CRS を持っていたので、簡単に答えることができたのである。

一般的に使われているデフォルトは Universal Transverse Mercator(UTM)で、地球を縦60個のくさびと横20個の緯度に分割した CRS のセットである。 UTM CRS で使用されている横メルカトル図法はコンフォーマルであるが、UTM ゾーンの中心から離れるにつれて面積や距離の歪みがひどくなる。 そのため、GIS ソフトウェア Manifold のドキュメントでは、UTM ゾーンを使用するプロジェクトの縦断範囲を中心子午線から6度までに制限することを提案している(出典: manifold.net)。 そのため、UTM は比較的狭い範囲でのアングル保存を重視する場合にのみ使用することを勧める。

地球上のほとんどの場所には UTM コードがあり、例えば「60H」 は R が発明されたニュージーランド北部を指している。 UTM EPSG コードは、北半球では32601から32660まで、南半球では32701から32760まで順次表示される。

その仕組みを説明するために、 ここにある通り、地球上の任意の地点に関連する EPSG コードを計算する関数 lonlat2UTM() を作ってみよう。

lonlat2UTM = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if(lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

次のコマンドは、この機能を利用して、オークランドとロンドンの UTM ゾーンと関連する EPSG コードを特定する。

lonlat2UTM(c(174.7, -36.9))
#> [1] 32760
lonlat2UTM(st_coordinates(london))
#> [1] 32630

現在、適切な CRS を選択するためのツールも用意されており、これには crssuggest パッケージ が含まれている。 このパッケージのメイン関数である suggest_crs() は、地理的な CRS を持つ空間オブジェクトを受け取り、与えられた領域に使用可能な投影 CRS のリストを返す。29 もう一つの便利なツールは、選択した場所とタイプに基づいて CRS をリストアップするウェブページhttps://jjimenezshaw.github.io/crs-explorer/。 重要な注意点:これらのツールは多くの場面で役立つが、適用する前に推奨される CRS の特性を知っておく必要がある。

適切な CRS がすぐにわからない場合、CRS の選択は、その後の地図や分析において保存することが最も重要である特性によって決められるべきだろう。 すべての CRS は、等面積、等距離、コンフォーマル(形状はそのまま)、またはそれらの妥協点の組み合わせである(「投影座標参照系」参照)。 ローカルパラメータを持つカスタム CRS を対象地域に合わせて作成し、単一の CRS がすべてのタスクに適合しないプロジェクトでは、複数の CRS を使用することができる。 測地線計算」は、CRS が適切でない場合のフォールバックを提供することができる(proj.org/geodesic.html を参照)。 どの投影 CRS を使っても、数百キロメートルに及ぶジオメトリでは、結果が正確でない可能性がある。

カスタム CRS を決定する際には、以下を勧める。30

  • カスタムローカル投影(原点の緯度・経度を調査地域の中心に設定)のランバート方位角等値投影( LAEA)、これはすべての場所で等面積投影だが数千キロメートル以上では形状が歪んでしまう。
  • ある地点とローカル投影の中心点との直線距離を具体的に正確に表す方位角等距離( AEQD)投影図
  • 数千 km に及ぶ地域のランバート共形円錐(LCC)投影。円錐は、距離と面積の特性がセカント線間で妥当となるように設定されている。
  • 極域のステレオグラフ(STERE)投影。ただし、中心から数千キロメートル離れた面積や距離の計算に頼らないように注意すること。

地域のデータセットに特化した投影 CRS を自動的に選択する方法として、調査地域の中心点に対して方位角等距離(AEQD)投影を作成することが考えられる。 これは、データセットの中心点に基づくメートル単位でカスタム CRS(EPSGコードなし)を作成するものである。 他のデータセットが作成したカスタム CRS と互換性がなく、数百キロメートルに及ぶ広範なデータセットに使用すると、結果が正確でなくなる可能性がある。

このセクションで説明する原則は、ベクタデータセットとラスタデータセットに等しく適用される。 しかし、CRS 変換のフィーチャには、それぞれの地理データモデルに特有のものがある。 ベクタデータ変換の特殊性は Section 7.7 で、ラスタ変換の特殊性は Section 7.8 で説明する。 次に、最後のセクションでは、カスタムマッププロジェクション(Section 7.9)を作成する方法を紹介する。

7.7 ベクタジオメトリの再投影

Chapter 2 は、ベクタ幾何学がいかに点で構成されているか、そしていかに点が直線や多角形などのより複雑なオブジェクトの基礎を形成しているかを示した。 つまり、ベクタの再投影は、直線や多角形の頂点となるこれらの点の座標を変換することになる。

Section 7.5 は、2つのオブジェクト間の距離を計算するために、少なくとも1つの sf オブジェクトを異なる CRS を持つ同等のオブジェクトに変換する必要がある例を含んでいる。

london2 = st_transform(london_geo, "EPSG:27700")

london の変換版ができたので、sf 関数 st_transform() を使って、2つのロンドン表現間の距離を求めることができる。31 londonlondon2 が 2 km 以上も離れているのは意外かもしれない。32

st_distance(london2, london_proj)
#> Units: [m]
#>      [,1]
#> [1,] 2018

CRS の問い合わせと再投影のための関数を、cycle_hire_osm を参照して以下に示す。これは spDatasf オブジェクトで、ロンドンで自転車をレンタルできる「ドッキングステーション」を表しているものである。 sf Section 7.1 オブジェクトの CRS は、関数 を使って問い合わせることができる。 st_crs() 出力は、座標系に関する情報を含む複数行のテキストとして印刷される。

st_crs(cycle_hire_osm)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCS["WGS 84",
#>     DATUM["WGS_1984",
#>         SPHEROID["WGS 84",6378137,298.257223563,
#>             AUTHORITY["EPSG","7030"]],
#>         AUTHORITY["EPSG","6326"]],
#>     PRIMEM["Greenwich",0,
#>         AUTHORITY["EPSG","8901"]],
#>     UNIT["degree",0.0174532925199433,
#>         AUTHORITY["EPSG","9122"]],
#>     AUTHORITY["EPSG","4326"]]

Section 7.3 で見たように、主な CRS コンポーネントである User inputwkt は単一の実体として出力される。st_crs() の出力は、実際には、次のコードチャンクの出力に示すように、inputwkt という単一の文字列という2つの要素を持つクラス crs の名前付きリストになっている。

crs_lnd = st_crs(london_geo)
class(crs_lnd)
#> [1] "crs"
names(crs_lnd)
#> [1] "input" "wkt"

Nameproj4stringepsg を含む追加の要素を $ 演算子で取り出すことができる(詳しくは ?st_crs と GDAL website の CRS and tranformation tutorial を参照)。

crs_lnd$Name
#> [1] "WGS 84"
crs_lnd$proj4string
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
crs_lnd$epsg
#> [1] 4326

Section 7.2 で述べたように、crs_lnd オブジェクトの $wkt 要素に格納された WKT 表現は、究極の真理の源である。 これは、前のコードチャンクの出力が、オブジェクトとその CRS の固有の属性ではなく、PROJ によって提供される wkt 表現からのクエリであることを意味する。

オブジェクトの CRS が変換されると、CRS の wktUser Input の両方の要素が変更される。 以下のコードチャンクでは、CRS を投影した新しいバージョンの cycle_hire_osm を作成する(簡潔さのため、CRS 出力の最初の4行のみを表示す)。

cycle_hire_osm_projected = st_transform(cycle_hire_osm, "EPSG:27700")
st_crs(cycle_hire_osm_projected)
#> Coordinate Reference System:
#>   User input: EPSG:27700 
#>   wkt:
#> PROJCRS["OSGB36 / British National Grid",
#> ...

この結果、オブジェクトは EPSG コード27700の新しい CRS を持つことになる。 しかし、この EPSG コードや他のコードの詳細を調べるにはどうしたらよいだろうか? ネットで検索するのも一つの方法である。

crs_lnd_new = st_crs("EPSG:27700")
crs_lnd_new$Name
#> [1] "OSGB 1936 / British National Grid"
crs_lnd_new$proj4string
#> [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
crs_lnd_new$epsg
#> [1] 27700

その結果、EPSG コード27700が英国ナショナルグリッドを表しており、「 EPSG 27700」とネット検索すれば出てくるような結果である。

コンソールで空間オブジェクトを表示すると、自動的にその座標参照系が返される。 明示的にアクセスして変更するには、st_crs() 関数、例えば st_crs(cycle_hire_osm) を使用する。

7.8 ラスタジオメトリの再投影

前節で説明した投影の概念は、ラスタにも同様に適用できる。 しかし、ベクタとラスタの再投影には重要な違いがある。 ベクタオブジェクトの変換は、すべての頂点の座標を変更することになるが、ラスタデータには当てはまらない。 ラスタは同じ大きさの矩形セル(度やメートルなどの地図単位で表現される)で構成されているため、通常、ピクセルの座標を個別に変換することは不可能である。 ラスタの再投影では、多くの場合、元のオブジェクトとは異なる列と行の数で、新しいラスタオブジェクトを作成する。 その後、新しいピクセルを適切な値で「埋める」ことができるように、属性を再推定する必要がある。 つまり、ラスタの再投影は、ラスタ範囲を別の CRS にベクタ再投影(Section 7.7)、リサンプリングによる新しいピクセル値の計算(Section 5.3.4)という2つの別々の空間処理として考えることができる。 したがって、ラスタデータとベクタデータの両方を使用する場合は、ラスタの再投影を避け、ベクタの再投影を行う方が良い場合がほとんどである。

通常のラスターの再投影は、ワーピングとも呼ばれている。 さらに、似たような操作で「変形」と呼ばれるものが2つある。 すべての値をリサンプリングするのではなく、すべての値をそのままに、ラスタセルごとに新しい座標を再計算し、グリッドジオメトリを変更する。 例えば、入力されたラスタ(正方形の格子)を曲線的な格子に変換することができる。 変換操作は、stars パッケージを用いて R で行うことができる。

ラスタの再投影処理は、terra パッケージの project() を使用する。 前節で紹介した st_transform() 関数と同様、project() は地理的オブジェクト(この場合はラスタデータセット)と何らかの CRS 表現を第二引数として受け取る。 余談だが、第2引数には、異なる CRS を持つ既存のラスタオブジェクトを指定することもできる。

ラスタ変換の例として、カテゴリデータと連続データを使ったものを見てみよう。 土地被覆データは、通常、カテゴリ化された地図で表現される。 nlcd.tif ファイルは、以下のコードチャンクの出力に示すように、NAD83 / UTM ゾーン 12N CRS の National Land Cover Database 2011 から取得した米国ユタ州の小領域の情報を提供する(出力の最初の行のみ示す)。

cat_raster = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
crs(cat_raster)
#> PROJCRS["NAD83 / UTM zone 12N",
#> ...

この地域では、8つの土地被覆クラスが区別された(NLCD2011の土地被覆クラスの全リストは、 mrlc.govで見ることができる)。

unique(cat_raster)
#>       levels
#> 1      Water
#> 2  Developed
#> 3     Barren
#> 4     Forest
#> 5  Shrubland
#> 6 Herbaceous
#> 7 Cultivated
#> 8   Wetlands

カテゴリ別ラスタを再投影する場合、推定値は元のラスタと同じでなければならない。 これは、各新規セルの値を入力ラスタの最も近いセル(中心)の値に設定する最近傍法(near)を使って行うことができる。 例えば、ウェブマッピングに適した地理的 CRS である WGS84 に cat_raster を再投影している。 まず、この CRS の PROJ 定義を取得する。これは、例えば、 http://spatialreference.org のウェブページを使用して行うことができる。 最後のステップは、project() 関数でラスタを再投影することである。カテゴリデータの場合は、最近傍法(near)を使用する。

cat_raster_wgs84 = project(cat_raster, "EPSG:4326", method = "near")

新しいオブジェクトの多くのプロパティは、Table 7.1 に示すように、列と行の数(したがってセルの数)、解像度(メートルから度に変換)、範囲など、以前のオブジェクトとは異なる(新しいカテゴリが作成されたのではなく、NA 値が追加されたため、カテゴリ数が 8 から 9 に増加していることに注意してみよう — 土地被覆クラスは維持されている)。

TABLE 7.1: オリジナル(‘cat_raster’)と 投影(‘cat_raster_wgs84’)のカテゴリーラスタ データセットの主要な属性。
CRS nrow ncol ncell resolution unique_categories
NAD83 1359 1073 1458207 31.5275 8
WGS84 1246 1244 1550024 0.0003 9

数値ラスタ(numeric またはこの場合は integer の値)の再投影もほぼ同じ手順で行う。 これは、Shuttle Radar Topography Mission (SRTM)spDataLarge にある srtm.tif で実証されており、WGS84 CRS による海抜メートル(標高)の高さを表している。

con_raster = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
crs(con_raster)
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

これから、このデータセットを投影型 CRS に再投影するが、カテゴリデータに適した最近傍法ではない。 その代わりに、元のラスタの4つの最近接セルに基づいて出力セルの値を計算するバイリニア法を使用する。33 投影されたデータセットの値は、これら4つのセルの値の距離加重平均である。 入力セルが出力セルの中心に近いほど、その重みは大きくなる。 以下のコマンドは、WGS 84 / UTM zone 12N を表すテキストストリングを作成し、この CRS にラスタを bilinear 方式で再投影するものである。

con_raster_ea = project(con_raster, "EPSG:32612", method = "bilinear")
crs(con_raster_ea)
#> [1] "PROJCRS[\"WGS 84 / UTM zone 12N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 12N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-111,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).\"],\n        BBOX[0,-114,84,-108]],\n    ID[\"EPSG\",32612]]"

数値変数のラスタ再投影は、セル数、解像度、範囲などの値や空間特性にも小さな変化をもたらす。 これらの変化は、Table 7.2 で実証されている。34:

TABLE 7.2: オリジナル(‘con_raster’)と投影(‘con_raster_ea’) 連続ラスタデータセット の主要な属性。
CRS nrow ncol ncell resolution mean
WGS84 457 465 212505 0.0008 1843
UTM zone 12N 515 422 217330 83.5334 1842
もちろん、2次元の地球投影の限界は、ラスターデータと同様にベクターデータにも当てはまる。 3つの空間特性(距離、面積、方向)のうち2つを維持することしかできない。 したがって、どのプロジェクションを選択するかは、目の前の課題によって決まる。 例えば、密度(1グリッドセルあたりの点数や1グリッドセルあたりの住民数)に興味がある場合は、等面積投影を使用する( Chapter 14 も参照)。

7.9 カスタム地図投影法

EPSG:4326 のような AUTHORITY:CODE の識別子で捕捉される確立された CRS は、多くのアプリケーションに適している。 しかし、場合によっては、代替的な予測を使用したり、カスタム CRS を作成したりすることが望ましい。 Section 7.6 は、カスタム CRS を使用する理由を述べ、いくつかの可能なアプローチを提示した。 ここでは、これらのアイデアを R で応用する方法を紹介する。

一つは、既存の CRS の WKT の定義を取り、その要素の一部を変更し、新しい定義を再投影に使用する方法である。 これは、空間ベクタでは st_crs()$wktst_transform() で、空間ラスタでは crs()project() で行うことができる。次の例では、zion オブジェクトをカスタム方位角等距離(AEQD)CRS に変換する例を示している。

zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))

カスタム AEQD CRS を使用するには、データセットの中心点の座標を度数で知る必要がある(地理的CRS)。 私たちの場合、zion 領域の重心を計算し、WGS84 に変換することでこの情報を抽出することができる。

zion_centr = st_centroid(zion)
zion_centr_wgs84 = st_transform(zion_centr, "EPSG:4326")
st_as_text(st_geometry(zion_centr_wgs84))
#> [1] "POINT (-113 37.3)"

次に、新たに得られた値を用いて、以下に示す方位角等距離(AEQD)CRS の WKT 定義を更新することができる。 "Central_Meridian" は重心の経度、"Latitude_Of_Origin" は緯度であることに注意しておこう。

my_wkt = 'PROJCS["Custom_AEQD",
 GEOGCS["GCS_WGS_1984",
  DATUM["WGS_1984",
   SPHEROID["WGS_1984",6378137.0,298.257223563]],
  PRIMEM["Greenwich",0.0],
  UNIT["Degree",0.0174532925199433]],
 PROJECTION["Azimuthal_Equidistant"],
 PARAMETER["Central_Meridian",-113.0263],
 PARAMETER["Latitude_Of_Origin",37.29818],
 UNIT["Meter",1.0]]'

この方法の最後のステップは、元のオブジェクト(zion)を新しいカスタム CRS(zion_aeqd)に変換することである。

zion_aeqd = st_transform(zion, my_wkt)

カスタムプロジェクションは、例えば Projection Wizard のウェブアプリケーションを使って、対話的に行うことも可能である (Šavrič, Jenny, and Jenny 2016)。 このサイトでは、データの空間的範囲と歪みのプロパティを選択すると、可能な投影のリストが返される。 また、このリストには投影の WKT 定義が含まれており、コピーして再投影に使用することができる。 WKT 文字列を用いたカスタム CRS 定義の作成については、Open Geospatial Consortium (2019) を参照。

PROJ 文字列は、Section 7.2 で述べた、投影、特に大きな地理的領域をカバーする幾何学に固有の制限を受け入れて、カスタム投影を作成するために使用することもできる。 多くの投影法が開発され、PROJ 文字列の +proj= の要素で設定することができる。 PROJ website だけでも数十のプロジェクトが詳細に記述されている。

面積の関係を維持したまま世界を地図化する場合、Figure 7.3 に示される Mollweide 投影が一般的で、賢明な選択となる (Jenny et al. 2017)。 この射影を使用するには、st_transform 関数の proj-string 要素 "+proj=moll" を使って指定する必要がある。

world_mollweide = st_transform(world, crs = "+proj=moll")
世界のモルワイデ図法。

FIGURE 7.3: 世界のモルワイデ図法。

世界地図を作成する際、すべての空間特性(面積、方向、距離)に対して歪みを最小化することが望まれることが多い。 これを実現するための代表的なプロジェクションとして、Figure 7.4 に示される ヴィンケル図法(第 3 図法) (Winkel Tripel Projections) がある。35 結果は、以下のコマンドで作成された。

world_wintri = st_transform(world, crs = "+proj=wintri")
世界のヴィンケル第 3 図法。

FIGURE 7.4: 世界のヴィンケル第 3 図法。

さらに、proj-string パラメータはほとんどの CRS 定義で変更可能であり、例えば +lon_0+lat_0 パラメータを使用して投影の中心を調整することができる。 以下のコードは、ニューヨークの経度と緯度を中心としたランバート方位等面積投影に座標を変換するものである(Figure 7.5)。

world_laea2 = st_transform(world,
                           crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=-74 +lat_0=40")
ニューヨークを中心とした世界のランバート方位等面積投影図。

FIGURE 7.5: ニューヨークを中心とした世界のランバート方位等面積投影図。

CRS の変更に関する詳しい情報は、Using PROJ のドキュメントに記載されている。

7.10 演習

E1. Create a new object called nz_wgs by transforming nz object into the WGS84 CRS.

  • Create an object of class crs for both and use this to query their CRSs.
  • With reference to the bounding box of each object, what units does each CRS use?
  • Remove the CRS from nz_wgs and plot the result: what is wrong with this map of New Zealand and why?

E2. Transform the world dataset to the transverse Mercator projection ("+proj=tmerc") and plot the result. What has changed and why? Try to transform it back into WGS 84 and plot the new object. Why does the new object differ from the original one?

E3. Transform the continuous raster (con_raster) into NAD83 / UTM zone 12N using the nearest neighbor interpolation method. What has changed? How does it influence the results?

E4. Transform the categorical raster (cat_raster) into WGS 84 using the bilinear interpolation method. What has changed? How does it influence the results?