【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

shimada yasu
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換すると
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

Mizutani Takayuki
島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

この部分のncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換すると
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss


_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

shimada yasu
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:
島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss

clipbord_R_res19.png (24K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

Mizutani Takayuki
水谷です。

ポイントのプロット部分を以下のようにadd=Tで重ねると大丈夫でした。
(par(new=T)で重ねると上手くいかないのかもしれません。わかりません)

r <- raster(ext, ncol=19, nrow=19)dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

2017年10月15日 16:52 shimada yasu <[hidden email]>:
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:

島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss


_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

shimada yasu
水谷様

日曜日にも係わらずありがとうございます。

やってみました。

(1)ポイントは、当初ダミーで設定したとおりの場所に表示されまし(添付)

(2)ただ、点はラスタの中心でなく、四隅に位置する・・ので、ラスタがどの点を代表しているのかが分からない。。。

 私が欲しかったのは、下記のような「ラスタのど真ん中に点がある」図なのですが。。。(考え方が間違っているかな。。)

-------------
|              |
|      +      |
|              |
------------     


しまだ

2017年10月15日 17:18 Mizutani Takayuki <[hidden email]>:
水谷です。

ポイントのプロット部分を以下のようにadd=Tで重ねると大丈夫でした。
(par(new=T)で重ねると上手くいかないのかもしれません。わかりません)

r <- raster(ext, ncol=19, nrow=19)dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

2017年10月15日 16:52 shimada yasu <[hidden email]>:
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:

島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss

clipbord_R_res19a.png (22K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

Mizutani Takayuki
水谷です。

この場合のraster関数は、第1引数を四隅の範囲として、
その範囲をncol、nrowの数で分割した解像度のラスタを作成します。
rasterize関数は、raster関数で作成した空のラスタにポイントの値を
入れ込むというイメージになります。
以下のコードではいかがでしょうか?

e<-extent(c(120.5, 140.5, 20.5, 40.5))
r <- raster(e, ncol=20, nrow=20)
dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

水谷

2017年10月15日 17:56 shimada yasu <[hidden email]>:
水谷様

日曜日にも係わらずありがとうございます。

やってみました。

(1)ポイントは、当初ダミーで設定したとおりの場所に表示されまし(添付)

(2)ただ、点はラスタの中心でなく、四隅に位置する・・ので、ラスタがどの点を代表しているのかが分からない。。。

 私が欲しかったのは、下記のような「ラスタのど真ん中に点がある」図なのですが。。。(考え方が間違っているかな。。)

-------------
|              |
|      +      |
|              |
------------     


しまだ

2017年10月15日 17:18 Mizutani Takayuki <[hidden email]>:

水谷です。

ポイントのプロット部分を以下のようにadd=Tで重ねると大丈夫でした。
(par(new=T)で重ねると上手くいかないのかもしれません。わかりません)

r <- raster(ext, ncol=19, nrow=19)dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

2017年10月15日 16:52 shimada yasu <[hidden email]>:
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:

島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss


_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

shimada yasu
水谷様

ありがとうございます。外なので、のちほど確認します。
なるほど、そのように考えればよいのですね。勉強になりました。

しまだ


2017/10/15 19:01 "Mizutani Takayuki" <[hidden email]>:
水谷です。

この場合のraster関数は、第1引数を四隅の範囲として、
その範囲をncol、nrowの数で分割した解像度のラスタを作成します。
rasterize関数は、raster関数で作成した空のラスタにポイントの値を
入れ込むというイメージになります。
以下のコードではいかがでしょうか?

e<-extent(c(120.5, 140.5, 20.5, 40.5))
r <- raster(e, ncol=20, nrow=20)
dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

水谷

2017年10月15日 17:56 shimada yasu <[hidden email]>:
水谷様

日曜日にも係わらずありがとうございます。

やってみました。

(1)ポイントは、当初ダミーで設定したとおりの場所に表示されまし(添付)

(2)ただ、点はラスタの中心でなく、四隅に位置する・・ので、ラスタがどの点を代表しているのかが分からない。。。

 私が欲しかったのは、下記のような「ラスタのど真ん中に点がある」図なのですが。。。(考え方が間違っているかな。。)

-------------
|              |
|      +      |
|              |
------------     


しまだ

2017年10月15日 17:18 Mizutani Takayuki <[hidden email]>:

水谷です。

ポイントのプロット部分を以下のようにadd=Tで重ねると大丈夫でした。
(par(new=T)で重ねると上手くいかないのかもしれません。わかりません)

r <- raster(ext, ncol=19, nrow=19)dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

2017年10月15日 16:52 shimada yasu <[hidden email]>:
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:

島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss


_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
Reply | Threaded
Open this post in threaded view
|

Re: 【R:教えて下さい】ポイントをラスタ変換すると中心がずれる

shimada yasu
水谷様

確認しました。素晴らしい、欲しかった図です。

同時にいくつかチェックすべき点を教えて頂きました。
今後はその点を確認しながら、作図したいと思います。

ありがとうございました。本当に助かりました。



2017年10月15日 20:12 shimada yasu <[hidden email]>:
水谷様

ありがとうございます。外なので、のちほど確認します。
なるほど、そのように考えればよいのですね。勉強になりました。

しまだ


2017/10/15 19:01 "Mizutani Takayuki" <[hidden email]>:
水谷です。

この場合のraster関数は、第1引数を四隅の範囲として、
その範囲をncol、nrowの数で分割した解像度のラスタを作成します。
rasterize関数は、raster関数で作成した空のラスタにポイントの値を
入れ込むというイメージになります。
以下のコードではいかがでしょうか?

e<-extent(c(120.5, 140.5, 20.5, 40.5))
r <- raster(e, ncol=20, nrow=20)
dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

水谷

2017年10月15日 17:56 shimada yasu <[hidden email]>:
水谷様

日曜日にも係わらずありがとうございます。

やってみました。

(1)ポイントは、当初ダミーで設定したとおりの場所に表示されまし(添付)

(2)ただ、点はラスタの中心でなく、四隅に位置する・・ので、ラスタがどの点を代表しているのかが分からない。。。

 私が欲しかったのは、下記のような「ラスタのど真ん中に点がある」図なのですが。。。(考え方が間違っているかな。。)

-------------
|              |
|      +      |
|              |
------------     


しまだ

2017年10月15日 17:18 Mizutani Takayuki <[hidden email]>:

水谷です。

ポイントのプロット部分を以下のようにadd=Tで重ねると大丈夫でした。
(par(new=T)で重ねると上手くいかないのかもしれません。わかりません)

r <- raster(ext, ncol=19, nrow=19)dd.map <- rasterize(dd, r, field=dd$z)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
plot(dd,add=T)

2017年10月15日 16:52 shimada yasu <[hidden email]>:
水谷様

ありがとうございました。

> ​この部分のncolとnrowは19になると思います。
> print(r)とするとresolutionで解像度が確認できます。

早速、やってみました。
解像度は1となりましたが、plotしたところ、中心にはなりせんでした・・・

-----
r <- raster(ext, ncol=19, nrow=19)
-----
class       : RasterLayer 
dimensions  : 19, 19, 361  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 121, 140, 21, 40  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
-----

2017年10月15日 15:30 Mizutani Takayuki <[hidden email]>:

島田さん

水谷です。
r <- raster(ext, ncol=20, nrow=20)

​​
この部分の
ncolとnrowは19になると思います。
print(r)とするとresolutionで解像度が確認できます。
この場合、ポイントはラスタの中心ではなく、ラスタの隅となります。

--
水谷貴行

2017年10月15日 12:03 shimada yasu <[hidden email]>:
すみません、どなたかご教授下さい。

地図表示をしたく、Rで格子状のポイントデータをラスタ変換する
ラスタの中心にポイントが来ないのですが、、、これは

(1)以下のスクリプトが間違っている?
(2)そもそも経度の距離は緯度に依存するから、平均的なラスタ解像度とは、当然一致しない?
(3)そのほか・・・

ご助言を頂けると幸いです。

RとGISも中途半端な理解なので、根本的に間違っているかも知れませんが、
数日考えてもどうにも解決できず投稿した次第です。

----ここから----
library(maptools)
library(rgdal)
library(raster)

rm(list=ls(all=TRUE))

#ダミーデータ(データフレーム)作成
x <- rep(seq(121,140,1),10)
y <- numeric(0);for ( i in 21:40 ){ y <- c(y, rep(i,20))}
z <- c(1:(length(x)*length(y)))

dd <- data.frame(x=x, y=y, z=z)

#座標を指定
coordinates(dd) <- ~x+y
#projectionを指定
proj4string(dd) <- CRS("+init=epsg:4326")

#領域設定
ext <- extent(dd)
r <- raster(ext, ncol=20, nrow=20)

#ラスタ変換
dd.map <- rasterize(dd, r, field=dd$z)

#左下隅を拡大(ラスタの中心に+がない・・・)
plot(dd.map, xlim=c(120,125), ylim=c(21,25))
par(new=T)
plot(dd, xlim=c(120,125), ylim=c(21,25))
------ここまで-----

--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss




--
島田泰夫/一般財団法人日本気象協会/環境・エネルギー事業部/環境影響評価室/〒170-6055豊島区東池袋3-1-1サンシャイン60-55F
/TEL:03-5958-8160 FAX:03-5958-8157 /[hidden email]
/携帯:080-8018-1567

_______________________________________________
OSGeoJapan-discuss mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss