Hola, mi trozo de código es:
# Primero creamos el mapa de toda la instrumentación
rm(list=ls()) # remove all the variables from the workspace
setwd("C:/Users/Marcos/Documents/TFM/Situacion_pluviometros/")
library(ggmap)
library(rgdal)
library(ggplot2)
library(scales)
library(maptools)
# Convertimos coordenadas de UTM a geográficas:
# Leemos el archivo de texto donde están las coordenadas UTM de los Pluviómetros
coord <-
read.delim("C:/Users/Marcos/Documents/TFM/Situacion_pluviometros/pluviometrosVenero.txt",
header=TRUE)[,4:5]
# Añadimos coordenadas del TDR, los dos limnímetros y el radar meteorológico
(figura 3.1 TFM_Carlos)
coord.todo <- rbind(coord,c(356489,4474915), c(360744,4470212),
c(360040,4471782),
c(359088,4474084))
rownames(coord.todo)=make.names(c(rep("Pluviómetros",7),
"Radar", "TDR", rep("Limnímetros",2)),
unique=TRUE)
coord.todo.spdf <- coord.todo # ".spdf": Class for spatial
attributes that have spatial point locations
colnames(coord.todo.spdf) <- c("x","y")
coordinates(coord.todo.spdf) <- ~x+y
# To retrieve the CRS for a spatial object: proj4string(x)
proj4string(coord.todo.spdf) <- CRS("+proj=utm +zone=30
ellps=WGS84")
# Ahora pasamos de un CRS a otro (pasamos de coordenadas UTM a geográficas)
coord.geog <- as.data.frame(spTransform(coord.todo.spdf,
CRS("+proj=longlat +datum=WGS84")))
colnames(coord.geog) <- c("lon", "lat")
# Creamos un data.frame llamado "datos" para poder pintar bien la
leyenda del mapa
datos <- coord.todo
datos$lon <- coord.geog$lon
datos$lat <- coord.geog$lat
datos$INSTRUMENTACIÓN <-
c("Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro",
"Radar
meteorológico","TDR","Limnímetro","Limnímetro")
# Ahora con la función "get_map" consultamos en Google Maps el mapa
# que queremos crear de nuestra zona de estudio
mapa <- get_map(location=c(lon=-4.66, lat=40.39), zoom=13,
maptype="terrain",
source="google")
nombres.pluvio <- c("Arromoro", "Atalaya",
"Dehesa", "Gasolinera", "Collado Morales",
"Peña Parda", "Trampalones")
# Cercamos el área de la cuenca
area <- readShapePoly("venero_Project.shp") # Lee el shape con
el area de la cuenca
proj4string(area) # Describes data's current
coordinate reference system
# To change to correct projection:
area <- spTransform(area, CRS("+proj=longlat +datum=WGS84"))
proj4string(area) <- "+proj=utm +zone=30 +north +ellps=WGS84
+datum=WGS84 +units=m +no_defs"
png("mapa_zona_estudio.png", units="in", width=9,
height=8.5, res=300)
# Vamos a pintar los distintos instrumentos de medida sobre el mapa (junto con
éste) y añadimos leyenda
ggmap(mapa, extend='device', legend="left",
base_layer=ggplot(datos, aes(x=lon, y=lat))) +
geom_point(data=datos, aes(shape = INSTRUMENTACIÓN, color=INSTRUMENTACIÓN,
fill=INSTRUMENTACIÓN), size=7.5) +
geom_text(data=coord.geog[1:7,], aes(x=coord.geog$lon[1:7],
y=coord.geog$lat[1:7], label=nombres.pluvio, hjust=0.2, vjust=1), size=5.2) +
scale_shape_manual(values=c(21,24,22,23))
geom_polygon(aes(x=long, y=lat, group=id), data=area, color="white",
alpha=.4, size=.2)
dev.off()
El error que me da la consola de R es:
Error in spTransform(xSP, CRSobj, ...) :
No transformation possible from NA reference system
¿Alguien me puede decir lo que hay que hacer para que ya no me de este error?
Muchas gracias de antemano.
Un saludo,
Marcos
[[alternative HTML version deleted]]