座標参照系と投影法
st_crs() と st_transform() を使って、CRS 間を正しく変換します。
「座標参照系と投影法」はCoddyKit上の無料R Academyレッスンです。 これはレッスン2/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
CRSが重要な理由
座標参照系(CRS)は、ジオメトリ列の数値が何を意味するのかをRに伝えます。たとえば、緯度・経度の度数、投影座標系におけるメートル、またはそれ以外の単位です。異なるCRSのデータを再投影せずに混在させると、空間結合がエラーを表示せずに失敗したり、距離計算が誤ったりします。
library(sf)
# Two sf objects in the same CRS: operations work
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# sf will error or warn when CRS do not match during operations
# This is the CRS metadata attached to the object
crs <- st_crs(nc)
cat('EPSG:', crs$epsg, '\n')
cat('Input is geographic (lon/lat)?', st_is_longlat(nc), '\n')
cat('Units:', crs$units_gdal, '\n')st_crs: CRS情報の読み取り
st_crs(sf_obj)は、EPSGコード、WKT文字列、Proj4文字列、単位情報を含むcrsオブジェクトを返します。整数のEPSGコードを直接渡して、CRSの定義を取得することもできます。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# Inspect the CRS object
crs <- st_crs(nc)
cat('Class:', class(crs), '\n')
cat('EPSG:', crs$epsg, '\n')
cat('Proj4:', substr(crs$proj4string, 1, 60), '...\n')
# Lookup a CRS by EPSG code
wgs84 <- st_crs(4326)
web_mercator <- st_crs(3857)
cat('\nWGS84 name:', wgs84$Name, '\n')
cat('Web Mercator name:', web_mercator$Name, '\n')EPSG:4326 — WGS84
EPSG:4326は、GPS座標の標準であるWorld Geodetic System 1984です。座標は10進度で表され、経度(X)は−180から180、緯度(Y)は−90から90の範囲です。これは地理座標系であり、投影座標系ではありません。
library(sf)
# Create points in WGS84 (GPS coordinates)
cities <- data.frame(
city = c('London', 'Tokyo', 'Sydney', 'New York'),
lon = c(-0.118, 139.692, 151.209, -74.006),
lat = c(51.509, 35.689, -33.869, 40.713)
)
cities_sf <- st_as_sf(cities,
coords = c('lon', 'lat'),
crs = 4326
)
cat('CRS:', st_crs(cities_sf)$Name, '\n')
cat('Units:', st_crs(cities_sf)$units_gdal, '\n')
# Naive distance in degrees (not meaningful for real distances)
d <- st_distance(cities_sf[1, ], cities_sf[2, ])
cat('London-Tokyo distance (degrees):', round(as.numeric(d), 2), '\n')EPSG:3857 — Web Mercator
EPSG:3857は、Google Maps、OpenStreetMap、Leafletで使われているWeb Mercator投影です。座標はメートルで表され、本初子午線と赤道を中心とします。極付近では面積が歪みますが、Webでの可視化には適しています。
library(sf)
cities <- data.frame(
city = c('London', 'Tokyo'),
lon = c(-0.118, 139.692),
lat = c(51.509, 35.689)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)
# Transform to Web Mercator
cities_merc <- st_transform(cities_sf, 3857)
cat('WGS84 coords (degrees):\n')
print(st_coordinates(cities_sf))
cat('\nWeb Mercator coords (metres):\n')
print(st_coordinates(cities_merc))
# Distance is now in metres
d <- st_distance(cities_merc[1, ], cities_merc[2, ])
cat('\nLondon-Tokyo distance (km):', round(as.numeric(d) / 1000, 0), '\n')st_set_crs: CRSの割り当て
st_set_crs(sf_obj, crs)は、座標を変換せずにCRSを割り当てます。CRSのメタデータが欠落している場合や誤っている場合にのみ使用してください。座標を新しいCRSに再投影する場合は、代わりにst_transform()を使用します。
library(sf)
# Create geometry with NO CRS
pts <- st_sfc(
st_point(c(-74.006, 40.713)),
st_point(c(-87.629, 41.878))
)
cat('CRS before:', is.na(st_crs(pts)), '(NA = no CRS)\n')
# Assign WGS84 (no coordinate transformation)
pts_wgs84 <- st_set_crs(pts, 4326)
cat('CRS after assignment:', st_crs(pts_wgs84)$epsg, '\n')
# Compare with st_transform (DOES transform coordinates)
pts_sf <- st_sf(id = 1:2, geometry = pts_wgs84)
pts_merc <- st_transform(pts_sf, 3857)
cat('Original lon:', st_coordinates(pts_sf)[1, 'X'], '\n')
cat('Mercator X: ', st_coordinates(pts_merc)[1, 'X'], '\n')st_transform: 再投影
st_transform(sf_obj, crs)は、すべてのジオメトリを現在のCRSから対象のCRSへ再投影し、完全な測地変換を適用します。crs引数には、EPSGの整数、Proj4文字列、またはWKT文字列を指定できます。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
cat('Original CRS:', st_crs(nc)$epsg, '\n')
# Re-project to WGS84
nc_wgs84 <- st_transform(nc, 4326)
cat('Transformed CRS:', st_crs(nc_wgs84)$epsg, '\n')
# Re-project to a local projected CRS (NC State Plane)
nc_stateplane <- st_transform(nc, 32119) # NAD83 / North Carolina
cat('State Plane CRS:', st_crs(nc_stateplane)$Name, '\n')
# Area calculation is more accurate in a projected CRS
area_deg <- st_area(nc[1, ])
area_m2 <- st_area(nc_stateplane[1, ])
cat('Area (degrees^2):', format(area_deg, big.mark = ','), '\n')
cat('Area (m^2): ', format(area_m2, big.mark = ','), '\n')適切な投影法の選択
投影法ごとに、最適化される性質が異なります。たとえば、正積図法(Mollweideなど)、正角図法(Mercator)、正距図法などがあります。局所的な分析ではメートル単位の投影CRSを使用し、全世界の概観にはWGS84を使用します。WebマップにはWeb Mercatorを使用します。
library(sf)
# Common EPSG codes cheat sheet
epsg_table <- data.frame(
EPSG = c(4326, 3857, 4269, 32637, 27700),
Name = c(
'WGS84 (GPS, global)',
'Web Mercator (web maps)',
'NAD83 (North America geographic)',
'UTM zone 37N (Middle East)',
'British National Grid (UK)'
),
Units = c('degrees', 'metres', 'degrees', 'metres', 'metres')
)
print(epsg_table, row.names = FALSE)
# Look up any EPSG online: https://epsg.io/<code>
cat('\nTip: use st_crs(epsg_code)$Name to verify a code\n')Proj4文字列
EPSGコードが標準になる前は、CRSはProj4文字列で定義されていました。これはPROJライブラリに渡す、簡潔なパラメータ文字列です。現在でもst_transform()で使用できますが、明確さと将来の互換性のためにEPSGコードが推奨されます。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# WGS84 as Proj4 string
wgs84_proj4 <- '+proj=longlat +datum=WGS84 +no_defs'
nc_wgs84 <- st_transform(nc, crs = wgs84_proj4)
cat('Using Proj4 string:\n', wgs84_proj4, '\n')
cat('Result CRS EPSG:', st_crs(nc_wgs84)$epsg, '\n')
# Compare: equivalent to st_transform(nc, 4326)
cat('Same as EPSG 4326?',
identical(st_crs(nc_wgs84), st_crs(st_transform(nc, 4326))), '\n')CRSの同一性の確認
空間結合やオーバーレイ操作を行う前に、st_crs(a) == st_crs(b)を使って両方のオブジェクトが同じCRSを共有していることを確認してください。異なる場合は、処理を続ける前に一方をもう一方に合わせて変換します。
library(sf)
# Two objects in different CRS
nc_nad27 <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84 <- st_transform(nc_nad27, 4326)
cat('Same CRS?', st_crs(nc_nad27) == st_crs(nc_wgs84), '\n')
# Always check before join/overlay
align_crs <- function(a, b) {
if (st_crs(a) != st_crs(b)) {
message('CRS mismatch! Transforming b to match a.')
b <- st_transform(b, st_crs(a))
}
list(a = a, b = b)
}
aligned <- align_crs(nc_nad27, nc_wgs84)
cat('After alignment - same CRS?',
st_crs(aligned$a) == st_crs(aligned$b), '\n')CRSを使った距離の計算
地理座標系(度数)上で計算した距離には意味がありません。まず正距または正積の投影CRSに再投影してください。または、sfに組み込まれた測地線距離の計算を使用することもできます。この計算はPROJライブラリを通じてWGS84を正しく処理します。
library(sf)
# City points in WGS84
cities <- data.frame(
city = c('Paris', 'Berlin', 'Madrid', 'Rome'),
lon = c(2.350, 13.404, -3.702, 12.496),
lat = c(48.865, 52.520, 40.417, 41.902)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)
# sf computes geodesic distances on WGS84 automatically
dist_matrix <- st_distance(cities_sf)
rownames(dist_matrix) <- cities$city
colnames(dist_matrix) <- cities$city
# Convert to km
dist_km <- round(as.numeric(dist_matrix) / 1000)
mat_km <- matrix(dist_km, 4, 4, dimnames = list(cities$city, cities$city))
print(mat_km)CRSを使った一連のワークフロー
完全なワークフローは次のとおりです。生のGPSデータを読み込み、WGS84を割り当て、正確な面積・距離計算のためにローカルなメートル単位のCRSへ再投影し、分析を実行してから、Webマッピング用の出力としてWGS84へ再投影します。
library(sf)
# Load and assign CRS
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84 <- st_transform(nc, 4326)
# Re-project to NC State Plane (metres) for accurate areas
nc_proj <- st_transform(nc, 32119)
# Compute area per county in km2
nc_proj$area_km2 <- as.numeric(st_area(nc_proj)) / 1e6
# Top 5 largest counties
top5 <- nc_proj[order(-nc_proj$area_km2), c('NAME', 'area_km2')]
cat('5 Largest NC Counties (km2):\n')
print(head(as.data.frame(top5)[, c('NAME', 'area_km2')], 5))
# Back to WGS84 for output
nc_output <- st_transform(nc_proj, 4326)
cat('\nOutput CRS:', st_crs(nc_output)$epsg, '\n')クイックチェック
現在WGS84(EPSG:4326)にあるポリゴンフィーチャについて、正確な面積(平方メートル)を計算する必要があります。どうすればよいでしょうか。
まとめ: CRSと投影法
重要なポイント:
st_crs(sf_obj)はCRSメタデータを読み取ります。CRSの指定にはEPSGコードが推奨されます- EPSG:4326 = WGS84(GPS用の度数)、EPSG:3857 = Web Mercator(Webマップ用のメートル)
st_set_crs()は変換せずにCRSを割り当て、st_transform()は座標を再投影します- 空間操作の前には、必ず
st_crs(a) == st_crs(b)でCRSが同じか確認します - 正確な距離や面積を求めるには、まずメートル単位の投影CRSへ再投影します
- Proj4文字列も引き続き使用できますが、EPSGコードのほうが簡潔で読みやすくなります
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# CRS pipeline
cat('Original:', st_crs(nc)$epsg, '\n')
nc_wgs84 <- st_transform(nc, 4326)
nc_proj <- st_transform(nc_wgs84, 32119)
cat('WGS84:', st_crs(nc_wgs84)$epsg, '\n')
cat('Projected:', st_crs(nc_proj)$epsg, '\n')AI チューターと学ぶ R — 無料
ブラウザでリアルコードを書いて実行し、24/7 の AI チューターから瞬時にサポートを受け、ウェブまたはアプリで続きから学習できます。
- コース
- 43
- レッスン
- 159
よくある質問
「座標参照系と投影法」レッスンは無料ですか?
はい。「座標参照系と投影法」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し(組み込みコードエディタと24時間対応のAIチューター)、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「座標参照系と投影法」で何を学びますか?
st_crs() と st_transform() を使って、CRS 間を正しく変換します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか?
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン2/4です。
「座標参照系と投影法」レッスンにはどのくらい時間がかかりますか?
ほとんどのCoddyKitレッスンは約5~10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか?
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。