Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #create raster file for Italy in right projection
- no2_2007Europe<-raster("D:/AELH/no2_2007/no2_07n/w001000.adf")
- no2_2007Italy<-crop(no2_2007Europe,extent(4e+06,5e+06,0,3e+06))
- 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"
- no2_2007Italy<-projectRaster(no2_2007Italy, crs="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs")
- #Load different trip files
- 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"))
- #Translate into spatial element with correct CRS
- latlong="+init=epsg:4326"
- CRS="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
- coords<-cbind(Longitude=as.numeric(as.character(trip_output$long)),Latitude=as.numeric(as.character(trip_output$lat)))
- trip_points<-SpatialPointsDataFrame(coords,trip_output[,-(3:4)],proj4string=CRS(latlong))
- trip_points<-spTransform(trip_points,CRS(AirMapCRS))
- ##To plot points##
- plot(trip_points,pch='.',col='darkblue',add=T)
- 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