############

##CALCULATE RURAL STATIONS WITHIN 1000 KM

#All rural stations within 1000 km are
#used to calculate the adjustment, with a weight that decreases linearly to zero at distance 1000 km.


##FUNCTIONS
   ### circledist
	#calculates great circle distance  with earth radius 6378.137 km
	#from http://en.wikipedia.org/wiki/Great-circle_distance 2nd formula	
	
   circledist =function(loc, lat,long,R=6372.795) {
	N=length(lat)	
	if(N==0) circledist=NA else {
	  pi180=pi/180;  
   	  x= abs(long -loc[2])
	  if (N>1) delta= apply( cbind(x %%360,(360-x)%%360),1,min)  *pi180  else delta= min (x %%360,(360-x)%%360) *pi180  
	  loc=loc*pi180; lat=lat*pi180; long=long*pi180
	  theta= 2* asin( sqrt(  sin( (lat- loc[1])/2 )^2 + cos(lat)*cos(loc[1])* (sin(delta/2))^2       ))
	  circledist=R*theta
	  }
	circledist
     }		
	
	circledist(loc= c(36.12,-86.67),lat=c(33.94),long= -118.40)
	 	#2887.258   #reported 2887

#GHCN INFO AT GISS: THE FULL LIST
	load("d:/climate/data/station/giss/giss.info.tab") 
	stations=giss.info; stations[1:3,]; names(stations)[3]="name"
	pi180=pi/180;R= 6372.795 #earth's radius in km

	stat_dist=rep(list(NA),7364)
	temp_rur=(stations$urban=="R")

	for (i in 1:7364) {
	   stations$dist=NA
	   lat0=stations$lat[i];long0=stations$long[i]
	   temp1= (abs (stations$lat-lat0) <10) 
		#initial screening only by latitude: must be within 10 deg latitude to be within 1000 km
	   test_angle= 10/sin(abs(lat0)*pi180)
	   x=	abs(stations$long -long0)
	   delta= apply( cbind(x %%360,(360-x)%%360),1,min) 
  	   temp2=(delta<test_angle) 
		#second screen: within 1000 km by longitude
	   stat_dist[[i]]=stations[temp1&temp2&temp_rur,c("id","name","long","lat","start_raw","end_raw","start_adj","end_adj","dist","pop","urban","lights")]
		#matrix of stations meeting first two screen tests
   	   stat_dist[[i]]$dist=circledist(loc=c(lat0,long0),lat=stat_dist[[i]]$lat,long=stat_dist[[i]]$long)
		#calculate circle distance for all stations in first screen
	   temp=(stat_dist[[i]]$dist<=1000)&!is.na(stat_dist[[i]]$dist)
	   stat_dist[[i]]=stat_dist[[i]][temp,]	
		#restrict to stations within 1000 km
	   stat_dist[[i]]=stat_dist[[i]][order(stat_dist[[i]][,"dist"]),]
		#order by increasing distance
	  }

	#save(stat_dist,file="d:/climate/data/station/giss/stat_dist.tab")
		#uploaded to www.climateaudit.org/data/station/giss/stat_dist.tab


 	for (i in 1:7364) {
   stat_dist[[i]]=stat_dist[[i]][order(stat_dist[[i]][,"dist"]),]