BeginPackage["Wikicode`NationalClimaticDataCenter`"]
StationsNearest::usage = "StationsNearest[{lat,lon}] finds the weather \
station closest to the given geocoordinates."
ISDFetch::usage = "ISDFetch[{usaf,wban,___},year] imports the ISD file \
corresponding to the given station (as returned by functions like \
StationsNearest) and year."
Begin["`Private`"]
(*Format specification:ftp://ftp.ncdc.noaa.gov/pub/data/ish/ish-\
format-document.pdf*)(*Examples:Import[isdfile,"ISD"] \
Import[isdfile,"ISD","Properties""{"Date","Time","Latitude",\
"Longitude","WindType"}]*)
ISDImport[f_, options___] := Module[{s, res}, s = OpenRead[f];
res = ReadList[s, String];
Close[s];
ISDParse[res, options]];
(*Additional data is not parsed,just returned in string form*)
$ISDBaseProperties =
Thread[{"USAF", "WBAN", "Date", "Time", "Source", "Latitude",
"Longitude", "ReportType", "Elevation", "CallLetter",
"QualityControl", "WindAngle", "WindDirection", "WindType",
"WindSpeed", "WindSpeedQuality", "CeilingHeight",
"CeilingQuality", "CeilingDetermination", "CAVOK",
"VisibilityDistance", "VisibilityQuality",
"VisibilityVariability", "VisibilityVariabilityQuality",
"Temperature", "TemperatureQuality", "DewPoint",
"DewPointQuality", "Pressure", "PressureQuality",
"Additional"} -> {(*{1,
4},*){5, 10}, {11, 15}, {16, 23}, {24, 27}, {28, 28}, {29,
34}, {35, 41}, {42, 46}, {47, 51}, {52, 56}, {57, 60}, {61,
63}, {64, 64}, {65, 65}, {66, 69}, {70, 70}, {71, 75}, {76,
76}, {77, 77}, {78, 78}, {79, 84}, {85, 85}, {86, 86}, {87,
87}, {88, 92}, {93, 93}, {94, 98}, {99, 99}, {100, 104}, {105,
105}, {106, -1}}];
$ISDStationProperties = {"USAF" -> {1, 6}, "WBAN" -> {8, 12},
"Name" -> {14, 42}, "Country" -> {44, 48}, "State" -> {50, 51},
"CallSign" -> {53, 57}, "Latitude" -> {59, 64},
"Longitude" -> {66, 72}, "Elevation" -> {74, 82},
"StartDate" -> {84, 91}, "EndDate" -> {93, 100}};
Options[ISDParse] = {"Properties" :> $ISDBaseProperties[[All, 1]]};
ISDParse[records_List, opts : OptionsPattern[]] :=
Module[{header =
SortBy[Intersection[$ISDBaseProperties[[All, 1]],
OptionValue["Properties"]], First[# /. $ISDBaseProperties] &],
data, nameToPos, p, m = Missing["NotAvailable"]},
nameToPos =
Append[Thread[
header -> Range[Length[header]]], _String :> (## &[])];
data = StringTake[records, header /. $ISDBaseProperties];
p = {"Date"} /. nameToPos;
p != {} && (p = p[[1]];
data[[All, p]] =
ToExpression@
StringTake[data[[All, p]], {{1, 4}, {5, 6}, {7, 8}}]);
p = {"Time"} /. nameToPos;
p != {} && (p = p[[1]];
data[[All, p]] =
ToExpression@StringTake[data[[All, p]], {{1, 2}, {3, 4}}]);
p = {"Latitude", "Longitude"} /. nameToPos;
p != {} && (data[[All, p]] = ToExpression@data[[All, p]]/1000.);
p = {"Elevation"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]] /. {9999 -> m});
p = {"WindAngle"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]] /. {999 -> m});
p = {"Temperature", "DewPoint", "WindSpeed"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]]/10. /. {999.9 -> m});
p = {"CeilingHeight"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]] /. {99999 -> m});
p = {"CAVOK"} /. nameToPos;
p != {} && (data[[All, p]] =
data[[All, p]] /. {"N" -> False, "Y" -> True, "9" -> m});
p = {"VisibilityDistance"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]] /. {999999 -> m});
p = {"VisibilityVariability"} /. nameToPos;
p != {} && (data[[All, p]] =
data[[All, p]] /. {"N" -> False, "V" -> True, "9" -> m});
p = {"Pressure"} /. nameToPos;
p != {} && (data[[All, p]] =
ToExpression@data[[All, p]]/10. /. {9999.9 -> m});
Thread[header -> Transpose[data]]]
ImportExport`RegisterImport["ISD", ISDImport]
stations =
Module[{m =
Missing["NotAvailable"]}, {# /. "999999" -> m, #2 /.
"99999" -> m, StringTrim@#3 /. "" -> m,
StringSplit@#4 /. {} -> m, StringTrim@#5 /. "" -> m,
StringTrim@#6 /. "" -> m,
If[#7 == "-99999", m, ToExpression@#7/1000.],
If[#8 == "-999999", m, ToExpression@#8/1000.],
If[StringTrim@#9 == "-99999", m, ToExpression@#9/10.],
If[#10 == "NO DATA ", m, DateList[#10][[;; 3]]],
If[#11 == "NO DATA ", m, DateList[#11][[;; 3]]]} & @@
StringTake[#, $ISDStationProperties[[All, 2]]] & /@
Select[ReadList[
StringToStream[
Import["http://www1.ncdc.noaa.gov/pub/data/inventories/ISH-\
HISTORY.TXT"]], String], StringLength@# == 100 &]];
StationsNearest =
Nearest[#[[7 ;; 8]] -> # & /@
Cases[stations, {__, _Real, _Real, _, _, _}],
DistanceFunction -> GeoDistance];
ISDFetch[{usaf_, wban_, ___}, year_] :=
Module[{m = Missing@"NotAvailable"},
Import["ftp://ftp.ncdc.noaa.gov/pub/data/noaa/" <> ToString@year <>
"/" <> If[usaf === m, "999999", usaf] <> "-" <>
If[wban === m, "99999", wban] <> "-" <> ToString@year <> ".gz",
"ISD"]]
End[]
EndPackage[]
