Programmēšana

Kā veikt telpisko analīzi R ar sf

Kur jūs balsojat? Kas jūs esat likumdevēji? Kāds ir jūsu pasta indekss? Šiem jautājumiem ir kaut kas ģeotelpiski kopīgs: atbilde ietver noteikšanu, kurā daudzstūrī punkts ietilpst.

Šādi aprēķini bieži tiek veikti ar specializētu ĢIS programmatūru. Bet tos ir viegli izdarīt arī R. Jums ir nepieciešamas trīs lietas:

  1. Veids, kā ģeokodēt adreses, lai atrastu platumu un garumu;
  2. Formfaili, kas iezīmē pasta indeksa daudzstūra robežas; un
  3. Sf pakete.

Ģeokodēšanai es parasti izmantoju geocod.io API. Tas ir bez maksas 2500 meklēšanas dienā, un tam ir jauka R pakete, taču, lai to izmantotu, jums ir nepieciešama (bezmaksas) API atslēga. Lai apietu šo raksta sarežģītību, es izmantošu bezmaksas atvērtā koda Open Street Map Nominatim API. Tam nav nepieciešama atslēga. Paketei tmaptools ir funkcija, geocode_OSM (), lai izmantotu šo API.

Ģeotelpisko datu importēšana un sagatavošana

Es izmantošu pakotnes sf, tmaptools, tmap un dplyr. Ja vēlaties sekot līdzi, ielādējiet katru ar pacman :: p_load () vai instalējiet nevienu, kas vēl nav jūsu sistēmā ar install.packages (), pēc tam ielādējiet katru ar bibliotēka ().

Šajā piemērā es izveidošu vektoru ar divām adresēm: mūsu biroju Framingemā (Masačūsetsā) un RStudio biroju Bostonā.

adreses <- c ("492 Old Connecticut Path, Framingham, MA",

"250 Northern Ave., Bostona, MA")

Ar geocode_OSM ģeokodēšana ir vienkārša. Rezultātus var redzēt, izdrukājot pirmās trīs kolonnas, ieskaitot platumu un garumu:

geocoded_addresses <- geocode_OSM (adreses)

drukāt (ģeokodētas_adreses [, 1: 3])

vaicājums lat lon

# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105

# 22 250 Northern Ave., Bostona, MA 42.34806 -71.03673

Ir vairāki veidi, kā iegūt pasta indeksa formas failus. Vienkāršākais, iespējams, ir ASV Tautas skaitīšanas biroja pasta indeksu tabulu apgabali, kas ir līdzīgi, ja ne tieši tādi paši kā ASV pasta pakalpojumu robežām.

ZCTA failu varat lejupielādēt tieši no ASV Tautas skaitīšanas biroja, taču tas ir fails visai valstij. Dariet to tikai tad, ja jums nav iebildumu par lielu datu failu.

Viena vieta, kur lejupielādēt ZCTA failu atsevišķai valstij, ir Census Reporter. Meklējiet visus datus pēc štatiem, piemēram, iedzīvotāju skaita, pēc tam pievienojiet pasta indeksu ģeogrāfijai un izvēlieties lejupielādes datus kā formas failu.

Es varētu manuāli lejupielādēt manu lejupielādēto failu, bet tas ir vieglāk R. Šeit es izmantoju bāzes R atvienot () funkciju lejupielādētajā failā un izpakojiet to projekta apakšdirektorijā ar nosaukumu ma_zip_shapefile. Tas junkpats = TRUE arguments saka, ka es nevēlos, lai atvienotu, pievienojot citu apakšdirektoriju, pamatojoties uz zip faila nosaukumu.

unzip ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpats = TRUE,

pārrakstīt = TRUE)

Ģeotelpiskais imports un analīze ar sf

Tagad beidzot kāds ģeotelpiskais darbs. Es importēšu formas failu R, izmantojot sf's st_read () funkciju.

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Lasīšanas slānis `acs2017_5yr_B01003_86000US02648 'no datu avota funkcijas un 4 lauki # ģeometrijas tips: MULTIPOLYGON # izmērs: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_def

Skrienot, esmu iekļāvis konsoles atbildi st_read () jo tur tiek parādīta informācija: epsg. Tas saka kāda koordinātu atskaites sistēma tika izmantota faila izveidošanai. Šeit tas bija 4326. Pārāk neiedziļinoties nezālēs, epsg būtībā norādakāda sistēma tika izmantota, lai trīsdimensiju globusa apgabalus - Zemi - pārveidotu par divdimensiju koordinātām (platumu un garumu). Tas ir svarīgi, jo pastāv a daudz dažādu koordinātu atskaites sistēmu. Es vēlos, lai manos pasta indeksa daudzstūros un adrešu punktos tiktu izmantots viens un tas pats, tāpēc tie kārtojas pareizi.

Piezīme: Šis fails nejauši ietver daudzstūri visam Masačūsetsas štatam, kas man nav vajadzīgs. Tāpēc es filtrēšu to Masačūsetsas rindu

pasta indekss_geo <- dplyr :: filtrs (pasta indekss,

nosaukums! = "Masačūsetsa")

Shapefile kartēšana ar tmap

Daudzstūra datu kartēšana nav nepieciešama, taču tā ir jauka mana formas faila pārbaude, lai redzētu, vai ģeometrija ir tāda, kādu es sagaidu. Izmantojot tmap’s, ātri varat uzzināt sf objekta grafiku qtm () (saīsne ātrai motīvu kartei) funkcija.

qtm (pasta indekss_geo) +

tm_legend (show = FALSE)

Ekrāni, kurus nošāva Šarona Mačla,

Izskatās, ka man tiešām ir Masačūsetsas ģeometrija ar daudzstūriem, kas varētu būt pasta indeksi.

Tālāk es vēlos izmantot ģeokodētās adreses datus. Pašlaik tas ir vienkāršs datu rāmis, taču tas jāpārvērš par sf ģeotelpisku objektu ar pareizo koordinātu sistēmu.

Mēs to varam izdarīt ar sf's st_as_sf () funkciju. (Piezīme: sf pakotnes funkcijas, kas darbojas ar telpiskajiem datiem, sākas ar st_, kas nozīmē “telpiskais” un “laicīgais”.)

st_as_sf () prasa vairākus argumentus. Zemāk esošajā kodā pirmais arguments ir pārveidojamais objekts - manas ģeokodētās adreses. Otrais argumentu vektors funkcijai norāda, kurām kolonnām ir vērtības x (garums) un y (platums). Trešais nosaka koordinātu atskaites sistēmu uz 4326, tāpēc tā ir tāda pati kā maniem pasta indeksa daudzstūriem.

point_geo <- st_as_sf (ģeokodētas_adreses,

koordinātas = c (x = "lon", y = "lat"),

crs = 4326)

Ģeotelpiskā pievienojas sf

Tagad, kad esmu izveidojis divas datu kopas, ar sf's ir viegli aprēķināt katras adreses pasta indeksu st_join () funkciju. Sintakse:

st_join (point_sf_object, polygon_sf_object, join = join_type)

Šajā piemērā es vēlos palaist st_join () uz ģeokodētajiem punktiem vispirms un pasta indeksa daudzstūriem otrajā vietā. Tas ir tā sauktais kreisās pievienošanās formāts: Viss tiek iekļauti punkti pirmajos datos (ģeokodētās adreses), bet tikai punkti, kas atbilst otrajos (pasta indekss) datos. Visbeidzot, mans pievienošanās veids ir st_within, tā kā es vēlos, lai mačs būtu punkti.

my_results <- st_join (point_geo, pasta indekss_geo,

pievienoties = st_within)

Tieši tā! Ja es skatos savus rezultātus, izdrukājot vairākas svarīgākās kolonnas, jūs redzēsiet, ka katrai adresei ir pasta indekss (slejā “nosaukums”).

drukāt (mani_rezultāti [, c ("vaicājums", "nosaukums", "ģeometrija")])

# Vienkārša funkciju kolekcija ar 2 funkcijām un 2 laukiem # ģeometrijas tips: POINT # izmērs: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # vaicājuma nosaukuma ģeometrija # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Kartēšanas punkti un daudzstūri ar tmap

Ja vēlaties kartēt punktus un daudzstūrus, šeit ir viens veids, kā to izdarīt, izmantojot tmap:

tm_shape (pasta indekss_geo) +

tm_fill () +

tm_shape (mani_rezultāti) +

tm_bubbles (col = "sarkans", izmērs = 0,25)

Ekrānuzņēmums: Sharon Machlis,

Vēlaties vairāk R padomu? Dodies uz lapu “Darīt vairāk ar R”!

$config[zx-auto] not found$config[zx-overlay] not found