269 lines
10 KiB
R
269 lines
10 KiB
R
# NS41 Satellite Data
|
|
# Michael D Cayton
|
|
# December 22, 2018
|
|
# Last update 12-23-18
|
|
|
|
### sum of col V81 through 99
|
|
### 21 not 240 remove
|
|
# number of rows, sum
|
|
# midified Bessal function of the second kind.
|
|
|
|
#### Load the libraries ###
|
|
|
|
library(tidyverse)
|
|
library(RCurl)
|
|
library(lubridate)
|
|
|
|
# Set the url for the ftp site working directory.
|
|
# This is where the origianal data is stored.
|
|
url <- 'ftp://ftp.ngdc.noaa.gov/STP/space-weather/satellite-data/satellite-systems/gps/data/ns41/'
|
|
|
|
# Create a vector of file names from the directory
|
|
# This will allow us to extract the dates and will be used to create a master file.
|
|
filename <- getURL(url, dirlistonly = TRUE)
|
|
filename
|
|
|
|
# Clean up the file list
|
|
destnames <- strsplit(filename, "\r*\n")[[1]]
|
|
destnames
|
|
|
|
# Create data frame for the data.
|
|
ns41 = NULL
|
|
ns41 <- tbl_df(ns41)
|
|
|
|
# For loop to access each file, combine into one dataframe, and add a column for the filename date
|
|
l <- length(destnames)
|
|
data <- NULL
|
|
i <- NULL
|
|
|
|
for (i in 1:l) {
|
|
data <- read.table(file = paste(url, destnames[i], sep=""))
|
|
data$Date <- as.Date(str_sub(destnames[i], start=6, end=11), "%y%m%d")
|
|
ns41 <- rbind(ns41, data)
|
|
print(paste(i, "of", l, "complete"))
|
|
}
|
|
|
|
# Save table.
|
|
write_csv(ns41, path="/Users/mdcay/Documents/Happy Dad/ns41.csv")
|
|
|
|
#Import raw data as needed.
|
|
ns41 <- read_csv(file="/Users/mdcay/Documents/Happy Dad/ns41.csv")
|
|
ns41 <- tbl_df(ns41)
|
|
|
|
# Select for V21 = 240
|
|
ns41 <- ns41 %>% filter(V21 == 240)
|
|
# Create table, grouped by unique week, and sum of V81 through V99
|
|
ns41_sums <- ns41 %>% group_by(Date) %>%
|
|
select(Date,V25, V82:V100) %>%
|
|
summarise(n=n(),
|
|
V24=sum(V25),
|
|
V81=sum(V82),
|
|
V82=sum(V83),
|
|
V83=sum(V84),
|
|
V84=sum(V85),
|
|
V85=sum(V86),
|
|
V86=sum(V87),
|
|
V87=sum(V88),
|
|
V88=sum(V89),
|
|
V89=sum(V90),
|
|
V90=sum(V91),
|
|
V91=sum(V92),
|
|
V92=sum(V93),
|
|
V93=sum(V94),
|
|
V94=sum(V95),
|
|
V95=sum(V96),
|
|
V96=sum(V97),
|
|
V97=sum(V98),
|
|
V98=sum(V99),
|
|
V99=sum(V100)
|
|
)
|
|
|
|
write_csv(ns41_sums, path='ns41_sums.csv')
|
|
ns41_sums <- read_csv(file="/Users/mcayton/Documents/Happy Dad/ns41_sums.csv")
|
|
# Create a list of the energy levels. These can be found in V100 through V130.
|
|
energies <- c(.6, .7, .8, .9, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 10)
|
|
|
|
# Use Modified Bessel function of the second kind. Each V81-V99 row is an interval and go with the energy list.
|
|
x <- ns41_sums %>% select(-Date, -n)
|
|
nu <- as.numeric(x[1,])
|
|
bessel <- besselK(nu, energies)
|
|
|
|
bessel
|
|
|
|
ggplot() +
|
|
geom_line(aes(nu, energies))
|
|
|
|
### sum of col V81 through 99
|
|
### 21 not 240 remove
|
|
# number of rows, sum
|
|
# midified Bessal function of the second kind. ### sum of col V81 through 99
|
|
### 21 not 240 remove
|
|
# number of rows, sum
|
|
|
|
### JSON with table info ###
|
|
# test <- fromJSON('{
|
|
# "decimal_day": { "DESCRIPTION":"Decimal Day of year",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "days",
|
|
# "START_COLUMN": 0
|
|
# },
|
|
#
|
|
# "Geographic_Latitude": { "DESCRIPTION":"Geographic_Latitude",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "degrees",
|
|
# "START_COLUMN": 1
|
|
# },
|
|
# "Geographic_Longitude": { "DESCRIPTION":"Geographic_Longitude",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "degrees",
|
|
# "START_COLUMN": 2
|
|
# },
|
|
# "Rad_Re": { "DESCRIPTION":"Radius (earth radii)",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "R_E",
|
|
# "START_COLUMN": 3
|
|
# },
|
|
# "rate_electron_measured": { "DESCRIPTION":"rate_electron_measured",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 4
|
|
# },
|
|
# "rate_proton_measured": { "DESCRIPTION":"rate_proton_measured",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 12
|
|
# },
|
|
# "collection_interval": { "DESCRIPTION":"collection_interval",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "seconds",
|
|
# "START_COLUMN": 20
|
|
# },
|
|
# "year": { "DESCRIPTION":"year eg 2015",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "years",
|
|
# "START_COLUMN": 21
|
|
# },
|
|
# "decimal_year": { "DESCRIPTION":"decimal_year contains fractional part of year",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "years",
|
|
# "START_COLUMN": 22
|
|
# },
|
|
# "svn_number": { "DESCRIPTION":"SVN number",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "none",
|
|
# "START_COLUMN": 23
|
|
# },
|
|
# "dropped_data": { "DESCRIPTION":"dropped_data=1 means ignore data",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "none",
|
|
# "START_COLUMN": 24
|
|
# },
|
|
# "b_coord_radius": { "DESCRIPTION":"radius from earths dipole axis",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "R_E",
|
|
# "START_COLUMN": 25
|
|
# },
|
|
# "b_coord_height": { "DESCRIPTION":"height above earths equatorial plane",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "R_E",
|
|
# "START_COLUMN": 26
|
|
# },
|
|
# "magnetic_longitude": { "DESCRIPTION":"magnetic longitude",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "degrees",
|
|
# "START_COLUMN": 27
|
|
# },
|
|
# "L_shell": { "DESCRIPTION":"L shell, dipole file, T89",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "R_E",
|
|
# "START_COLUMN": 28
|
|
# },
|
|
# "bfield_ratio": { "DESCRIPTION":"Bsatellite/Bequator",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "none",
|
|
# "START_COLUMN": 29
|
|
# },
|
|
# "local_time": { "DESCRIPTION":"Magnetic local time",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "h",
|
|
# "START_COLUMN": 30
|
|
# },
|
|
# "b_sattelite": { "DESCRIPTION":" field at satellite",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "gauss",
|
|
# "START_COLUMN": 31
|
|
# },
|
|
# "b_equator": { "DESCRIPTION":"B field at equator",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "gauss",
|
|
# "START_COLUMN": 32
|
|
# },
|
|
# "diffp": { "DESCRIPTION":"not sure what this is",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "gauss",
|
|
# "START_COLUMN": 33
|
|
# },
|
|
# "sigmap": { "DESCRIPTION":"not sure what this is",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "gauss",
|
|
# "START_COLUMN": 34
|
|
# },
|
|
# "electron_background": { "DESCRIPTION":"electron_background estimated",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 35
|
|
# },
|
|
# "proton_background": { "DESCRIPTION":"proton_background estimated",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 43
|
|
# },
|
|
# "proton_activity": { "DESCRIPTION":"proton_activity=1 if proton activity is occurring",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "none",
|
|
# "START_COLUMN": 51
|
|
# },
|
|
# "electron_temperature_fit": { "DESCRIPTION":"electron temperature from Maxwellian fit",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "MeV",
|
|
# "START_COLUMN": 52
|
|
# },
|
|
# "electron_density_fit": { "DESCRIPTION":"electron number density fit",
|
|
# "DIMENSION": [1],
|
|
# "UNITS": "cm^-3",
|
|
# "START_COLUMN": 53
|
|
# },
|
|
# "model_counts_electron_fit": { "DESCRIPTION":"E1-E8 rates from Maxwellian fit",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 54
|
|
# },
|
|
# "dtc_counts_electron": { "DESCRIPTION":"E1-E8 rates dead time corrected",
|
|
# "DIMENSION": [8],
|
|
# "UNITS": "hertz",
|
|
# "START_COLUMN": 62
|
|
# },
|
|
# "integral_flux_instrument": { "DESCRIPTION":"integral of electron flux fit above integral_flux_energy[i]",
|
|
# "DIMENSION": [30],
|
|
# "UNITS": "cm^-2sec^-1sr^-1",
|
|
# "START_COLUMN": 70
|
|
# },
|
|
# "integral_flux_energy": { "DESCRIPTION":"energies for integral_flux_instrument",
|
|
# "DIMENSION": [30],
|
|
# "UNITS": "MeV",
|
|
# "START_COLUMN": 100
|
|
# },
|
|
# "electron_diff_flux_energy": { "DESCRIPTION":"energies for the fluxes in electron_diff_flux_energy",
|
|
# "DIMENSION": [15],
|
|
# "UNITS": "MeV",
|
|
# "START_COLUMN": 130
|
|
# },
|
|
# "electron_diff_flux": { "DESCRIPTION":"electron flux at energies electron_diff_flux[i]",
|
|
# "DIMENSION": [15],
|
|
# "UNITS": "cm^-2sec^-1sr^-1MeV^-1",
|
|
# "START_COLUMN": 145
|
|
# }
|
|
# }')
|
|
|
|
|