Dipu
Dipu

Reputation: 123

Download, Plot map, and extract data in R

I downloaded a monthly data from [NASA data][1] and saved in .txt and .asc format. I am trying to plot and extract the data from the ASCII file, but unfortunately I am unable to do so. I tried the following:

1.

infile <- "OMI/L3feb09.txt"
data <- as.matrix(read.table(infile, skip = 3, header = FALSE, sep = "\t")) 
data[data == -9999] = NA
rr <- raster(data, crs = "+init=epsg:4326")
extent(rr) = c(179.375, 179.375+1.25*288, -59.5, -59.5+1*120)

Tried to extract for australia

adm <- getData("GADM", country="AUS", level=1)
rr = mask(rr, adm)
plot(rr)
library(rgdal)
r = raster("OMI/L3feb09.txt")
plot(r)
library(raster)
r = raster("OMI/L3feb09.txt")
plot(r)

4.Also tried,

df1 <- read.table("OMI/L3feb09.txt", skip = 11, header = FALSE, sep = "\t")

Tried the following from Stackoverflow link 1 Stackoverflow link 2

The problem is there are strings in the file in between number, such as "lat = -55.5"

Appreciate any kind of help. Thank you

[2]: https://stackoverflow.com/questions/42064943/opening-an-ascii-file-using-r

Upvotes: 0

Views: 343

Answers (1)

ahmathelte
ahmathelte

Reputation: 623

So, I downloaded one file and played around with it! It is not the best solution, however, I hope it can give you an idea.

library(stringr)
# read data
data<-read.csv("L3_tropo_ozone_column_oct04",header = FALSE, skip = 3,sep = "")

# this "" will seperate lat = -59.5  to 3 rows, and will be easier to remove.
#Also each row in the data frame constrained by 2 rows of "lat", represents #data on the later "lat".
lat_index<-which(data[,1]=="lat")

#you need the last row that contains data not "lat string 
lat_index<-lat_index-1

#define an empty array for results. 
result<-array(NA, dim = c(120,288),dimnames = list(lat=seq(-59.5,59.5,1),
                                                   lon=seq(-179.375,179.375,1.25)))

I assumed data -on 3 three digits- on each latituide is dividable by 3 resulting in 288, which equals the lon grid number. Correct me if I'm wrong.

# function to split a string into a vector in which each string has three letter/numbers 
split_n_parts<-function(input_string,n){
  
  # dislove it to many elements or by number 
  input_string_1<-unlist(str_extract_all(input_string,boundary("character")))
  
  output_string<-vector(length = length(input_string_1)/n)
  
  for ( x in 1:length(output_string)){
    
    output_string[x]<-paste0(input_string_1[c(x*3-2)],
                             input_string_1[c(x*3-1)],
                             input_string_1[c(x*3)])
  }
  
    return(as.numeric(output_string))
}

Here, the code loops, collects, write each lat data in the result array

# loop over  rows constrainted by 2 lats, process it and  assign to an array 
for (i in 1:length(lat_index)){
  
  if(i ==1){
    
    for(j in 1:lat_index[i]){
      
      if(j==1){
        
        row_j<-paste0(data[j,])
        
      }else{
        
        row_j<-paste0(row_j,data[j,])
      }
    }
    
  }else{
    ii<-i-1
    
    lower_limit<-lat_index[ii]+4
    upper_limit<-lat_index[i]
    
    for(j in lower_limit:upper_limit){
      
      if(j==lower_limit){
        
        row_j<-paste0(data[j,])
        
      }else{
        
        row_j<-paste0(row_j,data[j,])
      }
    }
    
  }

  result[i,]<-split_n_parts(row_j,3)
}

Here, is the final array and image

#plot as image 
image(result)

image!

EDIT: To continue the solution and put the end-result:

# because data is IN DOBSON UNITS X 10

result<-result/10

#melt to datafrome 
library(plyr)
result_df<-adply(result, c(1,2))

result_df$lat<-as.numeric(as.character(result_df$lat))
result_df$lon<-as.numeric(as.character(result_df$lon))
# plotting
library(maps)
library(ggplot2)
library(tidyverse)

world_map <- map_data("world")

#colors      
jet.colors <-colorRampPalette(c("white", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) 

ggplot() + 
  
  geom_raster(data=result_df,aes(fill=V1,x=lon,y=lat))+
  
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group),
               fill=NA, colour = "black")+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0))+
  
  scale_fill_gradientn(colors = jet.colors(7))

end-result

Upvotes: 1

Related Questions