地図XMLの公共座標系を持つ筆を数える②：地番（筆）の数を市区町村別に数える（Rを利用）

# はじめに
## 今回の狙い
前回の記事で、地図XMLの公共座標系割合（ファイル数ベース）を市区町村別に出しました。今回は、地図XMLの公共座標割合を、地番数ベースで計算したいと思います。その作業メモです。

↓前回の記事
https://qiita.com/T-ubu/items/391ad0197c23a22ebdc7

## 私の環境
Windows PCを使い、解析にはRを使いました。

* Windows 11 home
* R 4.3.2

## 対象
2023年度の地図XML

# 手順
## Step 1: 方法を考える
各地図XMLファイルに含まれる筆の数を数えるにはどうすればよいか考えます。
地図XMLファイルをGeoJSON等に変換してポリゴンの数を数えれば筆の数はわかります。

ありがたいことにG空間情報センターからダウンロードしてきた地図XMLのパッケージには、各ダウンロードパッケージごとのリストが格納されています。地番ごとにリストされていますし、筆界未定地であっても各地番がカウントされているので、ここのリストの行数は地番の件数（筆の件数）とみなしてもよいかと思います。**下のリストの1列目のファイル名ごとに行数を数えればよいということになります。**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/2bf4d366-3255-3928-7c45-ab7ddaefc379.png)

厳密にいうとV1、V2などで表現される「めがね地番（複数の土地が 1 筆で扱われる地番）」はリスト上で１つではないので、筆数と一致しません。W1等で表現される二重地番、複数の図郭に分割記載の分属管理地番、そして無地番のものもあるので本リストと地番数の差異は多少あるかと思いますが、地図XMLのポリゴンとしての筆数はこちらを数えればいいと思います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/a8d39afb-52ff-ddde-e9b8-69b8c4d7410d.png)

## Step 2: リストを準備する
2023年度の地図XMLダウンロードパッケージは2,005あります。それぞれのパッケージに含まれる2,005個のファイルを作業フォルダに集めます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/5936c8f3-764c-bfc2-436f-cab41f555a20.png)

## Step 3: 各ファイルの行数を数える
リストは2005の「市区町村-法務局」ごとに準備されており、その2005のパッケージごとに複数の地図XMLが含まれています。１つの地図XMLファイル（zipされているときは市区町村コード-法務省コード-通し番号.zipのファイル名）に複数の筆があるので、XMLファイルごと（最小単位のZIPファイルごと）に筆の数を数えます。

targetとして各ZIPパッケージのsearch-listを読み込み、その一覧からZIPファイル名（XMLファイル名）ごとに数を数えます（stat1のこと）。

```R
#ライブラリのインストール
library("dplyr") #データベースの操作をするためのパッケージ（filterなど）
library("stringr") #文字列処理のパッケージ

#対象ファイルリストの読み込み
listDataPath <- 'list'
dataList <- list.files(listDataPath, pattern="search-list.csv")

#各リストをターゲットとして、統計をだす処理をする。
for (i in 1:length(dataList)){
  if(exists("target")){rm(target)}
  target <- dataList[i]

  if(exists("shikuchoson")){rm(shikuchoson)}
  shikuchoson <- read.csv(paste(listDataPath,"/",target,sep=""),fileEncoding="cp932",header=TRUE,colClasses=rep("character",6))

  if(exists("stat1")){rm(stat1)}
  stat1 <- 
    shikuchoson %>% 
    group_by(ZIPファイル名) %>%
    summarise(n=n())
  print(stat1, n=nrow(stat1))

  for (j in 1:nrow(stat1)){
    write(paste(stat1[j,1],",",stat1[j,2],sep=""),"result0218.csv",append=TRUE)
  }
}
```

こんな感じで、各ZIPファイルに含まれる筆の数が出力されていきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/186d5e02-0dff-870c-a9e2-0cb4614d4911.png)


## Step 4: 各ファイルの座標系情報を統合
Step3までで、各XMLファイル中の筆数が分かったので、先日作ったファイル名と座標系の情報（下のもの）を結合させます。これでXMLファイル名、座標系、筆数のテーブルができるはずです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/6de9f504-c3fe-337e-49e5-5db2228d7560.png)
↑このファイルの作り方はこちら→ https://qiita.com/T-ubu/items/3a279d84d2bee134a9c7

統合相手になるStep３で作ったファイルはファイル名がZIPなので調整します。

```R
fudeKazu <- read.csv("result0218.csv",fileEncoding="UTF-8",header=FALSE,colClasses=c("character","numeric"))
zahyo <- read.csv("xml-projection-list.csv",fileEncoding="UTF-8",header=FALSE,colClasses=rep("character",2))

fudeKazu[,1] <- gsub(".zip",".xml",fudeKazu[,1])
result <- left_join(zahyo,fudeKazu,by=c("V1"="V1"))
write.csv(result, file = "./XMLprojectionNumber.csv")
```

これで、XMLファイル名、その座標系、それに含まれる筆数が出力できました。315,034のXMLファイルについての情報が出てきました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/4712dc3a-6a09-4e4e-93c5-184a0389a06f.png)

## Step 5: 市区町村ごとに座標系の割合を計算
今、市区町村コード-法務局コード-通番（例：01101-4300-1）とある番号ですが、これを最初の5桁の市区町村ごとに抜き出して座標系の割合を計算したいと思います。

＜流れのイメージ＞
1. 最初の5桁でレコードを抽出
↓
1. 座標系別に数を集計
↓
1. 割合を出す

まず、ラベル名が分かりにくかったので、少しタイトル名を変えます。ファイル名はMLprojectionNumber2.csv という名前にしておきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/1e1938f3-24ab-648b-4045-2a219674f925.png)

そうしてから以下のスクリプトをRで実行します。dplyrのライブラリを使います。
```R
#ライブラリ読み込み
library(dplyr)
#ファイル読み込み
motherTable <- read.csv("XMLprojectionNumber2.csv",fileEncoding="UTF-8",header=TRUE,colClasses=c("character","character","character","numeric"))
#余分なスペースを削除
motherTable[,3] <- gsub(" ", "", motherTable[,3] )

#市区町村コードを抜き出す、重複を除いてリストする
codes <- unique(str_sub(motherTable[,2],1,5)) #length(codes) は1992

#市区町村ごとのforループ
for( i in 1:length(codes)){
  if(exists("code")){rm(code)}
  code <- codes[i]
  # 対象市区町村のレコードだけを抽出
  if(exists("target")){rm(target)}
  target <- filter(motherTable,str_detect(file, paste(code,"-",sep="")))

  #統計をとる（変数を初期化してから）
  if(exists("stat1")){rm(stat1)}
  if(exists("all")){rm(all)}
  if(exists("nini")){rm(nini)}
  if(exists("kokyoRate")){rm(kokyoRate)}
  stat1 <- 
    target %>% 
    group_by(proj) %>%
    tally(count)
  #計算
  all <- sum(target[,4]) #すべての筆数
  if (nrow(stat1[stat1$proj=="任意座標系",1]) == 0){
    nini <- 0 
  } else {
    nini <- stat1[stat1$proj=="任意座標系",2] #任意座標系の筆数
  }
  kokyoRate <- round((all-nini)/all*100,2) #公共座標の割合（％）

  #出力
  write(paste(code,all,all-nini,kokyoRate,sep=","),"KokyoRate-shikuchoson.csv",append=TRUE) #市区町村コード、筆数、公共座標の筆数、割合（％）
}
```

こうして以下のように、**市区町村番号、筆数、公共座標の筆数、公共座標の割合**が出力されていきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/a0e11b0f-5fa7-1c31-5f0c-1fcaf80b4f3d.png)

## Step 6: 市区町村コードから市区町村名も追加
前回記事のStep 3と同様に市区町村コードから市区町村名を追加

```R
if(exists("result")){rm(result)}
if(exists("result2")){rm(result2)}
if(exists("code")){rm(code)}
result <- read.csv("KokyoRate-shikuchoson.csv",fileEncoding="UTF-8",header=FALSE,colClasses=c("character","numeric","numeric","numeric"))
code <- read.csv("code.csv",fileEncoding="UTF-8",header=FALSE,colClasses=rep("character",2))
result2 <- left_join(result,code,by=c("V1"="V1"))
write.csv(result2, file = "./KokyoRate-shikuchoson_with_names.csv")
```

こんな感じになりました
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/c9ceb818-6ea2-2df9-7213-9839449f6f87.png)



# まとめ
今回、市区町村ごとに地図XMLの筆数と、そのうちの公共座標の割合を数えることができました。
成果はGithubにも置いておきます。
https://github.com/ubukawa/chizu-xml1/blob/main/KokyoRate-shikuchoson_with_names.csv

都道府県別にも数えられると思います。地図XMLの地番区域ごとも方法を考えれば数えられます。

# 参考
地図XML仕様： https://www.moj.go.jp/content/000116464.pdf
前の記事：　https://qiita.com/T-ubu/items/391ad0197c23a22ebdc7