방명록

  1. 한창훈 2017.10.13 00:09 신고  수정/삭제  댓글쓰기

    선생님, 오래전부터 이 블로그에서 많은 도움받고 지금도 참조하고 있습니다.
    감사드립니다 광고 배너 눌고 갈께요 !!!

    • R Friend R_Friend 2017.10.13 00:18 신고  수정/삭제

      안녕하세요 한창훈님,
      블로그가 도움이 되었다니 기쁘네요.
      배너 광고 클릭해주셔서 감사합니다.
      주말에 커피숍가서 블로그 글 쓸 때 용돈으로 유용하게 사용하겠습니다. ^^

  2. 이상규 2017.10.11 11:02 신고  수정/삭제  댓글쓰기

    알프렌드님 안녕하세요! 연휴는 잘 보내셨나요 ㅎㅎ!

    제가 10억건정도의 레코드를 bind_row를 통해서 합치려고 합니다

    밑에처럼 코드를 작성했습니다.

    rprp <- data.frame()
    ## PM1 / 조건별 관찰자
    for(qq in 1:4036){
    for(ww in 3:12){
    for(ee in 1:3){
    attach(ko.in)
    aa <- filter(ko.in, investigator_id == qq, standard_num == ww, order_num == ee)
    detach(ko.in)


    ## 모니터링 데이터
    popo <- data.frame()


    for(l in 1:72){
    if(aa[l, 8] != 0){
    kk <- filter(PM1, food_name == aa[l,6], hour_num == aa[l, 7])
    zoro <- aa[l, 8] / 10
    for(oo in 1:zoro){
    popo <- rbind(popo, kk)
    }
    }
    }


    epep <- data.frame(popo$result_numeric)
    #bind로 합치기
    rprp <- bind_rows(rprp, epep)
    }
    }
    }


    근데 Error in if (aa[l, 8] != 0) { : missing value where TRUE/FALSE needed라는 오류가 떳는데
    위에 조건에 해당하는 데이터가 없는 것 같은데 해결방법이 궁금합니다

    • 이상규 2017.10.11 11:48 신고  수정/삭제

      oder_num 변수가 2일때
      standard_num가 3이 존재하지 않습니다.

      그래서 돌아가는 중에 aa에 데이터가 생기지 않는데 이런 데이터가 없을 때 무시하고 진행하는 방법이 궁금합니다 ㅜ!

    • R Friend R_Friend 2017.10.11 12:00 신고  수정/삭제

      tryCatch() 한번 구글링 해보실래요?

    • 이상규 2017.10.11 13:09 신고  수정/삭제

      넵! 공부해보고 오겠습니다

  3. Eureka 2017.10.07 22:18 신고  수정/삭제  댓글쓰기

    안녕하세요 운영자님! 추석은 잘 보내셨는지요!
    저는 아직 많은 진전을 이루지 못했습니다 ㅠㅠ 그저 구동이 되는 코드문을 찾아헤메는 것 같아요.
    이전에 여쭤보았던 코드문은 한줄한줄 풀어갈수록 더 복잡한 것만 같아, 다른 코드를 찾아서 구동해보고 있습니다..문제는 시간이 많지 않아서 큰일입니다. 다만, 제가 아래의 코드를 찾았고, 일단 구동은 되는데, 약간의 변형이 필요할 것 같지만 방법을 알 수 없습니다.

    우선 제가 가지고 있는 자료가 1948년 부터 2017년까지의 월 강수량을 포함하고 있는 "netCDF"형식의 파일인데, 이를 raster로 만들어 stack을 하기 위해

    data <- raster("prate.sfc.mon.mean.nc", varname="prate")
    data1 <- stack("prate.sfc.mon.mean.nc")

    위 두줄을 제가 임의로 추가하였습니다.

    궁극적으로 저는, 1,3,6,9,12 개월의 SPI (표준강우지수)를 산출하고자 하는데, 아래의 코드문에서 이렇게 개월을 설정할 수 있는 line이 있는지 의문입니다. 아무것도 변형하지 않고, 코드를 그대로 구동했을때 산출되는 값조차도 valid한 range에 있지는 않습니다만 output이 나오기는 하는 것 같습니다.

    1,3,6,8,12 SPI란 각 숫자에 해당하는 강수량이 해당 분석기간의 평균(강수량)에 비해 얼마나 떨어져 있는지를 살펴보는 것입니다. 따라서, 12SPI의 값이 valid range에 해당하는 -2에서 2중, -2에 가까울수록 강수량이 평균에 비해 부족하고, 따라서 가뭄이 있었다는 말이 됩니다.

    제가 사용한 코드의 원문과, 웹사이트를 첨부드리니 검토를 부탁드립니다...ㅠㅠ
    이 코드를 RStudio창에 그대로 복사했을경우 48번째에 해당하는 라인에.
    #replace Inf with the value that sets sd(spi)=1 or mean = 0 (minimises ff)
    라고 되어있어서, 이 부분이 변형할 수 있는 부분인 것 같아, 49번째에서
    spi.t[spi.t==Inf] <- optimize(ff,lower=0,upper=100)$minimum

    Inf에 정수 3을 예시로 대입해서 output을 보았더니, 값 자체가 산출이 안된 상태로 export 되었습니다.

    운영자님...확인을 부탁드립니다 ㅠ__ㅠ
    다시 한번 감사드립니다..!!

    본 코드는, http://joewheatley.net/wp-content/uploads/2011/10/spi_functions.txt
    여기에서 참고하였습니다..

    =========================================================

    library(raster)
    library(sp)
    library(rgdal)
    library(ncdf4)
    setwd("C:/test/")

    data <- raster("prate.sfc.mon.mean.nc", varname="prate")
    data1 <- stack("prate.sfc.mon.mean.nc")


    getPrecOnTimescale <- function(precipitation,k){
    # precipitation is a vector of monthly precipitation values
    # returns monthly precipation averaged over current month and prior k-1 months
    Nt <- length(precipitation)
    prec.k <- as.vector(sapply(seq(from=1, to=Nt),function(t) {tm <- max(t-k+1,1); mean(as.vector(precipitation[tm:t]))}))
    return(prec.k)
    }


    getSPIfromPrec <- function(precipitation){

    #takes a vector of precipitation values
    #and returns a vector of spi values


    Nt <- length(precipitation)
    #include full years data only
    years <- ceiling(Nt/12)
    full.years <- trunc(Nt/12)
    Nt.f <- full.years*12
    monthsInCurrentYear <- Nt %% 12
    #monthly analysis
    spi.o <- array(NA,c(years,12))
    for (m in 1:12)
    {
    mdata <- precipitation[seq(from=m, to=Nt, by=12)]
    # empirical cdf's for each month
    fit.cdf <- ecdf(mdata)
    #locate each data point on cdf by applying fit.cdf function for each location
    #cdf probabilites
    # these will be unformly distributed on interval 1/years
    cdfs <- as.vector(sapply(mdata,fit.cdf))
    #invert normal
    spi.t <- qnorm(cdfs)
    spi.tp <- spi.t[ spi.t != Inf] #drop Inf
    ff <- function(x) (1-sd(c(x,spi.tp)))^2
    #replace Inf with the value that sets sd(spi)=1 or mean = 0 (minimises ff)
    spi.t[spi.t==Inf] <- optimize(ff,lower=0,upper=100)$minimum
    # ensure mean is zero. spi.t is normally distributed with mean zero and sd approx 1
    spi.t <- spi.t - mean(spi.t)
    ifelse( !(monthsInCurrentYear==0),ifelse (m <=monthsInCurrentYear, spi.o[,m]<-spi.t, spi.o[,m]<- c(spi.t,0) ), spi.o[,m]<-spi.t)
    }

    spi <- array(0,c(Nt))

    for ( t in 1:Nt){
    month <- (t-1)%%12 + 1
    year <- 1 + trunc((t-1)/12)
    spi[t] <- spi.o[year,month]

    }
    return(spi)
    }

    output3 <- calc(data1, getSPIfromPrec)

    writeRaster(output3, filename = "C:/test/test6.tif", format="GTiff")

    • R Friend R_Friend 2017.10.08 13:27 신고  수정/삭제

      안녕하세요.

      제가 오늘은 교회에 하루 종일 있다가 저녁에 7시쯤 집에 들어가요. 오늘 밤에 R코드 보고 의견 드릴께요.

      핸드폰으로 얼핏 봐서는 소스 코드 문제보다는 input 데이터 포맷이 소스 코드에서 사용한 것과 다르기때문 이거나, 아니면 마지막에 사용자 정의 함수 적용에서 문제(아마도 getPrecOnTimescale() 함수 미적용?)가 있어서 output이 원하는 형태로 안나오는게 아닌가 의심이 듭니다.

      저녁에 집에 가서 PC로 좀 보고 의견 드릴께요.

    • Eureka 2017.10.08 16:45 신고  수정/삭제

      운영자님 주일에도 답변주시고 정말 감사합니다 ㅠ__ㅠ !!

      말씀해 주신대로 제가 사용한 input이 다릅니다. 원래의 코드에 사용된 data와 제가 사용한 자료가 "netCDF"파일로 형식은 같고, 개발자가 해놓은 것과 같은 자료를 사용했을 경우 output이 산출은 됩니다.

      다만,
      1) 이렇게 동일한 자료를 산출한다고 할 경우에도, 바로 윗글에서 여쭤본 바와 같이 48번째의 "Inf"에서 1,3,6,9,12와 같이 제가 원하는 개월 수 만큼을 묶을수는 없는 것 같습니다. 이 부분을 변형하는 것이 아닌 것 같아요..!

      2) 그리고 코드의 16번째 줄에 보면,
      #precipitation is a vector of monthly precipitation values#
      #returns monthly precipitation averaged over current month and prior k-1 month# 이렇게 되어있는데, 제가 삽입한 자요는 vector가 아닌 raster 형태였습니다.

      제가 임의로,
      data <- raster("prate.sfc.mon.mean.nc", varname="prate")
      data1 <- stack("prate.sfc.mon.mean.nc")

      위의 두줄을 추가하여 구동했을때 output은 나왔고, 원래 output이 가져야 할 값의 범위도 valid 했습니다.

      3) 그래서, 두가지 자료가 모두 "netCDF" 형식인데 왜 하나는 읽히고 하나는 안읽히는지 궁금해져서, 모두 raster를 stack으로 쌓은 후 어떻게 입력되어 있는지 보기 위해 R에서 그냥 읽어봤습니다.

      [A 자료의 형식]

      class : RasterStack
      dimensions : 2000, 7200, 14400000, 440 (nrow, ncol, ncell, nlayers)
      resolution : 0.05, 0.05 (x, y)
      extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
      coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
      names : X1981.01.01, X1981.02.01, X1981.03.01, X1981.04.01, X1981.05.01, X1981.06.01, X1981.07.01, X1981.08.01, X1981.09.01, X1981.10.01, X1981.11.01, X1981.12.01, X1982.01.01, X1982.02.01, X1982.03.01, ...

      [B 자료의 형식]
      class : RasterStack
      dimensions : 94, 192, 18048, 837 (nrow, ncol, ncell, nlayers)
      resolution : 1.875, 1.904129 (x, y)
      extent : -0.9375, 359.0625, -89.49406, 89.49406 (xmin, xmax, ymin, ymax)
      coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
      names : X1948.01.01.00.32.08, X1948.02.01.00.32.08, X1948.03.01.00.32.08, X1948.04.01.00.32.08, X1948.05.01.00.32.08,

      ==> 두 자료의 차이점이라고 한다면,

      - Dimension
      - Resolution (공간해상도가 거의 100배 가까이 차이가 있는 자료입니다)
      - extent (둘다 global data인데, 다르게 나타나네요...?)
      - names: A는 (년/월/일)이라면, B는 (년/월/일/(뒤에 알수없는..시간 표시같습니다)

      이렇게 되어있습니다.

      저는 궁극적으로는 A자료를 사용하고 싶은데, 코드 자체가 B의 자료를 사용하고 있다보니 output을 display 할 수가 없네요 ㅠㅠ

      만약 코드를 변형할 수 없도록 되어있다면...거기까지는 아직 생각할 수가 없습니다 ㅠㅠㅋㅋ현재까지 너무 많은 시도를 해봤고 이게 그나마 가능성이 많은 것이거든요..!

      항상 감사합니다 운영자님! ^^
      좋은 주일 되세요~

    • R Friend R_Friend 2017.10.08 23:46 신고  수정/삭제

      Eureka님, R script 올려주신대로 저도 똑같이 따라해봤는데요, 제일 마지막에 함수 적용하는 부분만 원래 소스대로 하니 SPI 지수가 원하는 결과가 나오는거 같습니다.

      (1) 데이터는 구글로 'prate.sfc.mon.mean.nc' 검색해서 아래 사이트에서 다운로드 받아서 사용했습니다.
      데이터 다운로드 => https://www.esrl.noaa.gov/psd/repository/entry/show?entryid=c1f514a3-11d7-47d7-a5f2-bd2978bddd0d


      (2) 데이터 읽어들이고 함수에 적용하는 것은 raster() 함수를 사용했습니다. (stack() 사용 안 함)

      data <- raster("prate.sfc.mon.mean.nc", varname="prate")


      (3) 아래 script는 원래 소스 그대로 사용하고 손을 안댔습니다. (1, 3, 6, 9, 12월 indexing 해오는 것과는 상관이 없구요, Infinite 값이 나왔을 경우 평균이 '0', 표준편차가 '1'이 나오도록 값을 변경해주는 함수입니다.)

      #replace Inf with the value that sets sd(spi)=1 or mean = 0 (minimises ff)
      spi.t[spi.t==Inf] <- optimize(ff, lower=0, upper=100)$minimum


      (4) sapply() 로 월 평균 벡터 구하고, SPI 구하는 함수를 data 에 적용하니 계산이 되네요. (시간이 무척 오래 걸리네요)

      spi <- sapply(1:12, function(i) getSPIfromPrec(getPrecOnTimescale(data,i)))


      (5) 1, 3, 6, 9, 12월은 칼럼 indexing 하시면 될거 같습니다.

      # spi of 1, 3, 6, 9, 12 each month
      spi_m01 <- spi[, 1] # spi of month 1
      spi_m03 <- spi[, 3] # spi of month 3
      spi_m06 <- spi[, 6] # spi of month 6
      spi_m09 <- spi[, 9] # spi of month 9
      spi_m12 <- spi[,12] # spi of month 12



      전체 R script는 아래와 같습니다.

      #=== R scirpt for SPI calculation ===
      rm(list=ls()) # clear all
      options(scipen=100) # nemeric notation (not exponential notation)

      # install and importing packages
      install.packages(c("raster", "sp", "rgdal", "ncdf4"))

      library(raster)
      library(sp)
      library(rgdal)
      library(ncdf4)

      # prate.sfc.mon.mean.nc (57.64MB) dataset
      # download from: https://www.esrl.noaa.gov/psd/repository/entry/show?entryid=c1f514a3-11d7-47d7-a5f2-bd2978bddd0d

      # set directory
      setwd("D:/admin/Documents/R/SPI/")

      data <- raster("prate.sfc.mon.mean.nc", varname="prate")

      str(data)
      head(data)

      getPrecOnTimescale <- function(precipitation,k){
      # precipitation is a vector of monthly precipitation values
      # returns monthly precipation averaged over current month and prior k-1 months
      Nt <- length(precipitation)
      prec.k <- as.vector(sapply(seq(from=1, to=Nt),function(t) {tm <- max(t-k+1,1); mean(as.vector(precipitation[tm:t]))}))
      return(prec.k)
      }


      getSPIfromPrec <- function(precipitation){

      #takes a vector of precipitation values
      #and returns a vector of spi values
      Nt <- length(precipitation)

      #include full years data only
      years <- ceiling(Nt/12)
      full.years <- trunc(Nt/12)
      Nt.f <- full.years*12
      monthsInCurrentYear <- Nt %% 12

      #monthly analysis
      spi.o <- array(NA,c(years,12))
      for (m in 1:12)
      {
      mdata <- precipitation[seq(from=m, to=Nt, by=12)]
      # empirical cdf's for each month
      fit.cdf <- ecdf(mdata)

      # locate each data point on cdf by applying fit.cdf function for each location
      # cdf probabilites
      # these will be unformly distributed on interval 1/years
      cdfs <- as.vector(sapply(mdata, fit.cdf))

      #invert normal
      spi.t <- qnorm(cdfs)
      spi.tp <- spi.t[ spi.t != Inf] #drop Inf
      ff <- function(x) (1-sd(c(x, spi.tp)))^2

      #replace Inf with the value that sets sd(spi)=1 or mean = 0 (minimises ff)
      spi.t[spi.t==Inf] <- optimize(ff, lower=0, upper=100)$minimum

      # ensure mean is zero. spi.t is normally distributed with mean zero and sd approx 1
      spi.t <- spi.t - mean(spi.t)
      ifelse( !(monthsInCurrentYear==0),ifelse (m <=monthsInCurrentYear, spi.o[,m]<-spi.t, spi.o[,m]<- c(spi.t,0) ), spi.o[,m]<-spi.t)
      }

      spi <- array(0,c(Nt))

      for ( t in 1:Nt){
      month <- (t-1)%%12 + 1
      year <- 1 + trunc((t-1)/12)
      spi[t] <- spi.o[year,month]

      }
      return(spi)
      }

      spi <- sapply(1:12, function(i) getSPIfromPrec(getPrecOnTimescale(data,i)))

      str(spi) # Large matrix [1:18048, 1:12]
      head(spi)

      # spi of 1, 3, 6, 9, 12 each month
      spi_m01 <- spi[, 1] # spi of month 1
      spi_m03 <- spi[, 3] # spi of month 3
      spi_m06 <- spi[, 6] # spi of month 6
      spi_m09 <- spi[, 9] # spi of month 9
      spi_m12 <- spi[,12] # spi of month 12

    • Eureka 2017.10.09 09:58 신고  수정/삭제

      운영자님 정말 감사합니다 ^^ 이렇게 직접 시도까지 해주시구요!

      저도 원본 데이터를 사용했을시에 큰 문제는 없었습니다. 다만 문제라면, 공간해상도가 200km 정도 되는 자료여서, 우리나라가 정사각형 격자로 한 5개 정도로 표현이 된다는 것이죠ㅠㅠ너무 rough해서 의미를 찾기 힘들것 같아 자료의 사용이 어렵습니다.

      그래서 대안으로 찾은것이 약 10km짜리 netcdf 파일인데요, 안에 대체 어떻게 다르게 되어있는지 영 구동이 되지를 않습니다 ㅜㅜ 어제 글에 말씀드린대로 dimension 등의 부분에서 차이가 있고, 년/월/일 형식으로 표시된 것이 차이인데요. 혹시 시간이 되신다면 아래의 자료를 한번 봐주실 수 있으신지요.ㅜㅜ부탁드립니다! 용량이 커서 죄송해요 ㅜ 해상도가 10배 이상 좋은 자료이다보니, 위의
      데이터보다 기간이 짧아도 용량이 훨씬 더 큰 것 같습니다.
      매번 큰 도움과 친절한 설명에 감사드려요!! ㅠㅠ

      ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0/global_monthly/netcdf/

    • R Friend R_Friend 2017.10.10 00:38 신고  수정/삭제

      Eureka님,

      (1) 알려주신 경로 'ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0/global_monthly/netcdf/'에서 데이터를 다운로드하려고 하니 '5일' 남았다는 메시지가 뜹니다. 데이터 사이즈는 5.8G로 무척 큰데 서버 네트웍 속도가 너무 느린거 같습니다. 아무래도 제가 데이터를 다운로드 받아서 확인해보는 것은 불가능할 거 같습니다.
      ㅜ.ㅜ


      (2) 방명록에 남겨주신 현재 상황에 대한 설명만 봐서는 input 데이터 포맷 문제인거 같습니다.

      두번째 댓글에서 A자료, B자료 비교를 해주셨는데요, A/B 자료 모두 'class : RasterStack' 라고 되어 있는데요, 이는 stack("prate.sfc.mon.mean.nc") 함수를 사용해서 데이터를 불어왔을 때의 class 입니다.

      대신에 raster("prate.sfc.mon.mean.nc", varname="prate") 함수를 써서 데이터를 불러오면 'class : RasterLayer'가 됩니다.

      => stack() 함수 대신에 raster() 로 데이터를 불러와서 input 으로 사용해보시지요.

      => 제일 마지막에 사용자정의함수
      적용하실 때 output3 <- calc(data1, getSPIfromPrec) 이렇게 하셨었는데요(즉, getPrecOnTimescale() 함수 미적용),

      spi <- sapply(1:12, function(i) getSPIfromPrec(getPrecOnTimescale(data,i))) 처럼 'getPrecOnTimescale() 과 getSPIfromPrec() 을 모두 적용해보시지요.

      => 1, 3, 6, 9, 12월은 아래처럼 indexing 하면 될거 같습니다.
      spi_m01 <- spi[, 1] # spi of month 1
      spi_m03 <- spi[, 3] # spi of month 3

      위에 안내해드린 가설이 만약 원하시는 output을 만들어내지 못하면 제가 더이상 도와드릴 수가 없을 것 같습니다.

      부디 잘 해결되면 좋겠네요.

    • Eureka 2017.10.13 11:17 신고  수정/삭제

      운영자님 답변 감사드립니다. ^^

      아무래도 더 공부를 해보아야 할 것 같아요!
      안에 형식이 어떤 식으로 차이가 있을지는 모르겠지만, nc가 아닌 래스터로 읽어들였을 때에도 구동될 수 있도록 해보아야 겠어요..ㅠ

      매번 감사합니다. 너무 유용해서 이제는 정말 하루에 꼭 한번 습관처럼 들리게 되네요 ^^

      다른 포스팅 보고 또 많은 학습이 될 것 같습니다. 좋은 내용 감사합니다!

  4. 안효준 2017.10.05 22:48 신고  수정/삭제  댓글쓰기

    좋은 자료 정말 감사합니다!!!

  5. Eureka 2017.10.05 20:15 신고  수정/삭제  댓글쓰기

    안녕하세요
    멀게만 생각했던 R에 대해 너무 좋은 정보과 설명, 업데이트 해주셔서 정말 감사드립니다!

    • 2017.09.06 18:18  수정/삭제

      비밀댓글입니다

    • R Friend R_Friend 2017.09.06 18:21 신고  수정/삭제

      안녕하세요 Eureca님,
      R 사용하신다니 반갑습니다.

      비밀 댓글로 샘플 데이터 100여개 정도 남겨주시면 R script 한번 짜보겠습니다.

      제가 오늘 야근이라서 집에 밤에 퇴근 후에 답변드릴수 있습니다. 오늘 안될수도 있구요.

      사례금은 안주셔도 돼구요, 가끔씩 블로그에 광고 눌러주시면 커피값 용돈으로 잘 쓰겠습니다. ^^

    • 2017.09.06 19:16  수정/삭제

      비밀댓글입니다

    • R Friend R_Friend 2017.09.07 08:55 신고  수정/삭제

      어제 밤에 R코드 반절 정도 짰습니다. 오늘 밤에 나머지 짜서 올려놓겠습니다. 어제 야근하기도 했고, 며칠전 주문한 고성능 데스크탑이 새로 와서 그거 새로 세팅하다보니 밤이 너무 늦어서 다 못짰네요.

      오늘도 퇴근 후에 수업듣는게 있어서 집에 늦게 들어가요. 그래서 R코드 마저 짜는게 밤에 많이 늦을거 같으니 기다리지 마시구요, 내일 오전에 확인해보시면 좋겠습니다.

    • 2017.09.07 10:23  수정/삭제

      비밀댓글입니다

    • R Friend R_Friend 2017.09.08 00:53 신고  수정/삭제

      Eureka 님, 아래 R script 참고하세요. 한번 실행시켜 보시고, 혹시 에러나면 댓글 남겨주세요. 시각화는 일단 히트맵으러 했는데요, 혹시 지도에 좌표별 매쉬업해야하는거면 구글링해서 시각화 R code 찾아보시면 있을거예요. 논문 잘 마무리 하시길 바랄께요. ^^

      #======================================================
      # correlation analysis of monthly rainfall and vitality
      #======================================================

      # installing and loading packages
      install.packages(c("reshape", "dplyr", "ggplot2", "DBI"))
      library(reshape)
      library(dplyr)
      library(ggplot2)

      # Importing txt dataset
      rainfall <- read.table("D:/R/rainfall/rainfall.txt", sep = "\t", header = TRUE)
      vitality <- read.table("D:/R/rainfall/vitality.txt", sep = "\t", header = TRUE)


      #---- rainfall
      # Reshaping data structure
      rainfall_colname <- colnames(rainfall)

      rainfall_melt <- melt(data = rainfall,
      id.vars = c('XCoord', 'YCoord'),
      measure.vars = rainfall_colname[3:length(rainfall_colname)])

      # Changing the variable name
      rainfall_melt <- dplyr::rename(rainfall_melt,
      year = variable,
      rainfall = value)

      #---- vitality
      # Reshaping data structure
      vitality_colname <- colnames(vitality)

      vitality_melt <- melt(data = vitality,
      id.vars = c('XCoord', 'YCoord'),
      measure.vars = vitality_colname[3:length(vitality_colname)])

      # Changing the variable name
      vitality_melt <- dplyr::rename(vitality_melt,
      year = variable,
      vitality = value)




      # Merging 'rainfall_melt' and 'vitality_melt' : Inner Join
      rainfall_vitality <- merge(x = rainfall_melt, y = vitality_melt,
      by = c('XCoord', 'YCoord', 'year'))

      # Delete temp objects
      rm(rainfall_colname, rainfall_melt, vitality_colname, vitality_melt)


      # Making XCoord & YCoord key
      rainfall_vitality <- transform(rainfall_vitality,
      XY_Coord = paste(XCoord, YCoord, sep="_"))

      # Making XY_Coord list
      XY_Coord_list <- distinct(rainfall_vitality, XY_Coord)


      #------ function and loop
      rainfall_vitality_cor <- data.frame() # blank data.frame

      for(i in 1:nrow(XY_Coord_list)){
      # subset per XY_Coord
      rainfall_vitality_sub <- subset(rainfall_vitality,
      subset = (XY_Coord == XY_Coord_list[i, 1]))

      # pearson correlation coefficients
      rainfall_vitality_sub_cor <- cor(rainfall_vitality_sub[,c('rainfall', 'vitality')],
      use = "complete.obs",
      method = "pearson")

      # selecting XCoord, YCoord, corr_coef
      rainfall_vitality_sub_cor_2 <- cbind(rainfall_vitality_sub[1,c('XCoord', 'YCoord')],
      corr_coef = rainfall_vitality_sub_cor[1,2])

      # binding sub result to 'rainfall_vitality_cor' data.frame
      rainfall_vitality_cor <- rbind(rainfall_vitality_cor, rainfall_vitality_sub_cor_2)

      # delete temp objects
      rm(rainfall_vitality_sub, rainfall_vitality_sub_cor, rainfall_vitality_sub_cor_2)

      print(paste0(i, " / ", nrow(XY_Coord_list))) # check progresss
      }

      str(rainfall_vitality_cor)


      #----------- Visualization
      # Bubble chart using ggplot2
      ggplot(rainfall_vitality_cor, aes(x = XCoord, y = YCoord, fill=corr_coef)) +
      geom_tile() +
      ggtitle("Heat map of correlation b/w rainfall and vitality")

    • 2017.09.08 21:21  수정/삭제

      비밀댓글입니다

  6. 2017.09.27 12:58  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend R_Friend 2017.09.09 22:44 신고  수정/삭제

      안녕하세요 Eureka 님,
      이미 R 상당히 잘 사용하고 계시네요.

      중간에 멘트 적어놓으신 것처럼 'N-timescale' 객체가 없기 때문에 에러가 난거네요.

      spi.m <- as.list(rep(0,N_timescale))

      코드 앞에다가 'N_timescale' 에 값을 할당해주시면 됩니다.

    • Eureka 2017.09.26 18:27 신고  수정/삭제

      안녕하세요 운영자님 ^^
      지지난주부터 내내 다른일로 바쁘다가, 다시 원래 업무로 복귀하게 되었습니다 ㅠ__ㅠ

      답변에 감사드립니다!

      음, 저번에 질문드린 부분에서부터 잘 이해가 가지 않는데요...코드 앞에다가 'N_timescale'에 값을 할당하는 것은 무엇인지요...? 'N_timescale'자리에, 그냥 숫자를 입력하는 것인가요?

      spi.m <- as.list(rep(0,3))

      이런식으로, 말그대로 숫자 값을 입력하는 것인지요. 또 한번 여쭤봅니다! 항상, 감사드립니다.

    • Eureka 2017.09.26 18:53 신고  수정/삭제

      운영자님,
      N_timescale에 값을 대입하는 것이 맞는것 같아요!

      일단, 1,3,6,9,12 개월에 대한 누적치를 구하려고 하는데요,

      이 값을...마지막 plot 형태가 아닌,

      1) X,Y 좌표 값을 가진 csv 형태로도 추출할 수 있는 것인가요~?
      만약할 수 있다면, 코드 상에서

      "#plot" 이전에, csv 형식문을 넣어야 할 것 같은데 SPI 값은 어떤 것을 통해 추출될 수 있는것인지요~?


      2) 또한 plotting을 할 때, [colorRampPalette]를 사용하고 싶지 않습니다. 최종적으로 나오는 것이 시계열적인 SPI의 값 변이만을 보여주고 있는 것 같은데, 저는 이것을 "공간적"으로 보고자 하거든요~

      따라서 저는, 입력자료로 사용된 자료가 가지고 있는 X,Y 좌표를 그대로 가져와 [지도]의 형태로 plotting 하려고 합니다. 저는 매년/월에 해당하는 SPI값을 해당 X,Y 위치에 따라 산출하고 싶으며, 이 자료를 csv의 형태로 가져와 통계분석을 하려고 합니다. 이러한 방법이 가능한 것인지요?

      답변을 부탁드립니다. 항상, 감사드립니다. ^^

    • R Friend R_Friend 2017.09.30 21:23 신고  수정/삭제

      제 포스팅의 "R 프로그래밍" 카테고리에서 for loop 혹은 사용자 정의함수 글들 참고하세요.

      write.table() 함수 안에 파일 이름을 paste() 함수를 사용해서 년/월 구분자를 자동으로 부여하게 하면 되겠네요.

      지도에 매핑하는건 저도 안해봤는데요, 구글링 해보시면 나올거 같습니다.

      혹시 RShiny 할 줄 아시면 웹으로 애플리케이션으로 해서 조회할 때마다 interactive하게 그래프 바뀌게 하면 편할거 같습니다.

  7. 이상규 2017.09.20 10:33 신고  수정/삭제  댓글쓰기

    안녕하세요 ~ 알프렌드님.. 또 찾아왔습니다. 핳핳..

    약 1000개의 plot을 저장해야 하는데 문제가 생겨서 해결이 안되서 여쭤보러 왔습니다

    데이터는
    x name
    1 A B
    2 A B
    1 A B

    이런식으로 되어 있고 인덱싱을 통해 name에 A B 처럼 띄어쓰기가 되어있습니다.
    제가 작성한 코드는

    png(paste0(item[1, 2], ".png"))로 작성했는데

    띄어쓰기가 없는 name은 문제없이 작성 되는데 이건 어떻게 해결해야 할지 모르겠네요 ㅜㅜ
    도와주세요!

    • R Friend R_Friend 2017.09.20 11:34 신고  수정/삭제

      제가 문제를 잘 이해한건지 모르겠는데요, 띄여쓰기가 에러를 유발한다면 공백을 없애주면 되지 않을까요?

      df <- trnasform(df, name_2 = gsub(" ", "", name))

    • 이상규 2017.09.20 12:47 신고  수정/삭제

      DB랑 연동해서 사용하려다보니 공백도 그대로 있어야 해서요 ㅜㅠㅜ
      조언해주셔서 감사합니다

  8. 이상규 2017.09.15 14:54 신고  수정/삭제  댓글쓰기

    안녕하세요. 알프렌드님

    for문을 작성했고 문제없이 잘 돌아가는데
    for문 앞에 function을 입력해서 함수로 만들었는데 실행은 되긴 하는데 결과값이 안나오네요..

    무슨 문제일까요 ㅜㅜ?

  9. 이상규 2017.09.14 11:31 신고  수정/삭제  댓글쓰기

    안녕하세요. 여러 방법으로 시도 해봤지만 결국 못해서 질문드리러 왔습니다.

    데이터의 형태가

    #Option1 : 20, 30
    #Option2 : 15, 17
    -0.4000, 0.
    -0.1200, 0.
    -0.3300, 0.

    이런식으로 이뤄져 있고 모두 한 컬럼에 있습니다.

    '#'으로 시작되는 옵션행을 모두 삭제하고 나머지 부분만 남기고 싶은데 안되더라구요ㅜㅜ

    도움을 요청하고 싶습니다

    • 이상규 2017.09.14 11:55 신고  수정/삭제

      추가적으로 질문 있습니다.
      위와 똑같은 데이터에서

      kk[31,1] 이런식으로 인덱싱을 했는데
      [1] -0.3800, 0.
      4133 Leverls : -0.01000, 0. -0.202000, 0. -0.03000, 0, ..... 이런식으로 나오는데
      다른 복합적인 데이터형태인지.. 궁금합니다. 간략히 설명해주시면 열심히 공부하겠습니다.

      매번 많이 배우고 갑니다. 감사합니다

    • R Friend R_Friend 2017.09.14 12:35 신고  수정/삭제

      제가 오늘 회사일이 너무 바빠서요, 퇴근 후에 집에 가서 밤에 답변 달아드릴께요.

      두번째 질문은 데이터형이 factor여서 그렇습니다. as.numeric()으로 데이터형태 변경해보세요.

    • 이상규 2017.09.14 17:35 신고  수정/삭제

      감사합니다!! 잘 해결했습니다!!
      추가적인 질문으로
      lm함수로 회귀분석을 진행한 후
      abline으로 회귀선을 그리고
      회귀식을 text로 입력하고 싶습니다.

      수동적인 방법으로 강제로 입력하는 법은 알고있는데 혹시 여러 함수를 입력해야 할 때 간단하게 입력하는 방법이 있을까요?

    • R Friend R_Friend 2017.09.15 15:04 신고  수정/삭제

      회귀모형 적합시킨 모형 이름을 lm.fit 이라고 했을 때요,

      회귀계수 인덱싱은
      lm.coeffs <- coefficients(lm.fit)
      으로 하구요,

      회귀변수 인덱상은
      lm.vars <- names(lm.fit)
      이런긱으로 하면 돼요.

      위의 두 개의 객체(회귀계수, 변수 리스트)를 가지고 순서대로 인덱싱해서 회귀식을 입력하는 함수를 짜시면 됩니다.

  10. rstudy 2017.09.08 21:03 신고  수정/삭제  댓글쓰기

    안녕하세요,r friend 블로그의 도움을 많이 받고 있는 학생입니다.
    r 데이터 처리에 대해 궁금한 것이 있습니다.
    텍스트를 구분하기 위해
    strsplit을 이용해 a,b,c를 a b c로 나눴습니다. 이 과정에서 strsplit을 사용하기 위해 as.character를 사용했습니다. 그리고 split한 텍스트를 리스트 형태로 저장을 하려고 합니다.
    코드는 이렇습니다.

    for(i in 1:n){
    index <- as.character(DF[i, 12])
    split<-as.list(strsplit(genre_index, ",")
    }

    그런데 저는 데이터 분석을 위해 character형이 아닌 다른 자료형으로 다시 변환이 필요한 상황입니다. 이런 경우 어떤 코딩을 추가해줘야할지 궁금합니다..

    근본적인 질문을 드리자면, 저는 한 튜플 안에 한 가지 데이터만 들어간다는 데이터베이스의 전제조건에서 해결이 되지 않는 궁금함이 있어, 리스트 자료형에 주목했습니다.
    각 튜플 안에 list자료형으로 a,b,c를 각각 한 튜플에 넣어주면 데이터 분석을 할 때 여러 속성에 대한 개별적인 분석이 가능하다고 생각했습니다. 이러한 접근 방법이 맞는지 궁금합니다.
    예를 들어 음식 데이터베이스에
    1row: 김치볶음밥, 재료 column: 김치, 양파,마늘,햄, 계란
    2row: 김밥 재료 column: 김, 밥, 햄, 오이, 계란

    이런식의 데이터베이스에 햄과 계란이 중복돼있듯이, 각각의 재료들에 대한 데이터를 추출해서 분석하려고 합니다.
    재료 테이블을 따로 만들지 않고 list 자료형을 이용해 데이터 분석을 하는 것이 맞는지 궁금합니다.

    • R Friend R_Friend 2017.09.09 23:05 신고  수정/삭제

      안녕하세요.

      (1) DataFrame 의 문자열 변수를 strsplit 함수를 사용해서 구분하는 방법은

      data.frame(do.call('rbind', strsplit(as.character(df$var), split='delimeter', fixed=T)))

      을 사용해보세요. (참고 포스팅 ☞ http://rfriend.tistory.com/37 )


      (2) 분석 목적, 분석 기법/알고리즘이 무엇이냐에 따라서 데이터 구조가 달라질텐데요, 일반적인 supervised learning 모델을 사용할 것이라면 각 변수별로 1, 0 의 코드값을 가진 가변수(dummy variable)을 가진 DataFrame을 만들어서 사용합니다. 가령, '햄'이 포함되어 있으면 '1', 없으면 '0', '계란'이 포함되어 있으면 '1', 없으면 '0' 으로 해서 데이터 프레임을 만드는 것이지요.

      [데이터프레임 예시]
      ========================
      food kimchi ham egg cumumber
      김치볶음밥 1 1 1 0
      김밥 0 1 1 1
      ========================


      혹시 연관성 분석하려는 것이라면 sparse matrices 형태로 저장을 합니다.
      (transaction이 일어나 상품 (즉, 값이 '1'인 경우) 에 대해서만 저장해서 저장 공간 절약 위하여...)
      연관성 분석을 위해 데이터셋 구성하는 방법은 아래 블로그 포스팅 참고하세요.

      http://rfriend.tistory.com/193

    • 2017.09.10 16:26  수정/삭제

      비밀댓글입니다

    • R Friend R_Friend 2017.09.11 14:56 신고  수정/삭제

      말씀해주신 분석 목적이라면

      (1) 영화 장르별로 변수를 만들어서 각 영화마다 해당되면 1, 안되면 0으로 해서 가변수 만들고 나서
      ( rfrienf.tistory.com/57 참고)

      (2) 관객수를 목적변수로, 장르별 볌수를 설명변수로 해서 회귀모형 적합 시키면 되겠네요. Stepwise method로 젼수섬택 하시구요.