Advertisement
Guest User

Untitled

a guest
Jun 27th, 2017
56
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.10 KB | None | 0 0
  1. #create raster file for Italy in right projection
  2. no2_2007Europe<-raster("D:/AELH/no2_2007/no2_07n/w001000.adf")
  3. no2_2007Italy<-crop(no2_2007Europe,extent(4e+06,5e+06,0,3e+06))
  4. projection(no2_2007Italy)<-"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
  5. no2_2007Italy<-projectRaster(no2_2007Italy, crs="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs")
  6.  
  7. #Load different trip files
  8. trip_output<-read.csv("D:/AELH/Raw data/ROMA/1",header=F,sep=";",col.names=c("trip_id","time_stamp","lat","long","vehicle_type","speed"))
  9.  
  10.  
  11. #Translate into spatial element with correct CRS
  12. latlong="+init=epsg:4326"
  13. CRS="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
  14.  
  15. coords<-cbind(Longitude=as.numeric(as.character(trip_output$long)),Latitude=as.numeric(as.character(trip_output$lat)))
  16. trip_points<-SpatialPointsDataFrame(coords,trip_output[,-(3:4)],proj4string=CRS(latlong))
  17.  
  18. trip_points<-spTransform(trip_points,CRS(AirMapCRS))
  19.  
  20. ##To plot points##
  21. plot(trip_points,pch='.',col='darkblue',add=T)
  22. plot(no2_2007Italy,axes=F,border=rgb(0.8,0.8,0.8),add=T)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement