2

I did some spatial interpolations in R and I need a support for those.

Looking for shapefiles online the best I could find was this and it looks really horrible, the same area in ggplot looks like this. I spent a lot of time on this, to no avail, I found shapefiles with any kind of informations, but nothing that has a good definition of the shore.

EDIT my code

rm(list=ls())
setwd("C:/Users/irene/Desktop/FIN_adm")
graphics.off()
Finland <- readOGR(".","FIN_adm0")
#OGR data source with driver: ESRI Shapefile 
#Source: ".", layer: "FIN_adm0"
#with 1 features and 70 fields
#Feature type: wkbMultiPolygon with 2 dimensions

#prj file says GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

setwd("C:/tesi/Rtoshiba/dati/3 downtown")
data44 <- read.table("ch44.txt", header=T)
head(data44)
locations<-unique(cbind(data44[,1], data44[,2]))
#plot(locations)

media.stat44=vector(length=109)
media.rssi44=vector(length=109)

for (i in 1:109){
  temp<-data44[data44[,1]==locations[i,1] & data44[,2]==locations[i,2],]#location i
  media.stat44[i]<-mean(temp$stat)
  media.rssi44[i]<-mean(temp$rssi)
}

temp44<-cbind(locations, media.stat44)
xy <- temp44[,1:2]
df44 <- as.data.frame(temp44[,-1:-2])
SPDF44 <- SpatialPointsDataFrame(coords=xy, data=df44)


bbox = c(min(data44$lon), min(data44$lat),
         max(data44$lon), max(data44$lat))

predgrid <- expand.grid(lon=seq(from=bbox[1], to=bbox[3], length.out=400), 
                        lat=seq(from=bbox[2], to=bbox[4], length.out=400))

aa<-as.data.frame(rep(NA,160000))  
dimnames(aa)[[2]]<-"stat"
loc<-SpatialPointsDataFrame(coords=predgrid, data=aa)

ciccia44<-as.data.frame(df44)
dimnames(ciccia44)[[2]]<-"stat"
SPDFgri44<-SpatialPointsDataFrame(coords=locations, data=ciccia44)

stima44 <- idw(stat~1, SPDFgri44, loc, idp = 1.0, debug.level = 1)

library(raster)
rast44 <- raster(stima44)

f44<-as.data.frame(stima44)
# create spatial points data frame
spg <- f44
coordinates(spg) <- ~ x + y
# coerce to SpatialPixelsDataFrame
gridded(spg) <- TRUE
# coerce to raster
rasterDF <- raster(spg)
rasterDF
x11()
plot(rasterDF)

# setwd("C:/tesi/GIS_data/FIN_adm")
# #setwd("C:/Users/irene/Desktop/fi")
# 
# Finland <- readShapePoly("FIN_adm1.shp")
# Finland
# #str(Finland)
# unique(Finland@data$VARNAME_1)
# Finland <- Finland[Finland$VARNAME_1 == "Etelä-Suomen lääni|Södra Finlands län",]
vv<-fortify(Finland)

#prova subset sia giusto
#x11()
#plot(vv[,1], vv[,2])

projection(Finland) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")

elevation.sub <- crop(rasterDF, extent(Finland))

elevation.sub <- mask(elevation.sub, Finland)

x11()
plot(elevation.sub)
plot(Finland, add = TRUE)

I downloaded (between others) the shapefile suggested in the answer, selected the tile I am interested in getting this result Code I used:

Finland2 <- readOGR(".","helsinki")


class       : SpatialPolygonsDataFrame 
features    : 1331 
extent      : 2704914, 2905588, 8392151, 8516090  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : error, tile_x, tile_y 
min values  :     0,    227,    283 
max values  :     0,    228,    284 
unique(Finland2@data$tile_x)
unique(Finland2@data$tile_y)

I select the tile I am interested into

Finland3 <- Finland2 [Finland2$tile_x == 227 &  Finland2$tile_y == 284, ]

And project it

Finland3<-spTransform(Finland3, CRS("+proj=longlat +ellps=WGS84 +no_defs+ellps=WGS84 +towgs84=0,0,0 "))

x11()
plot(Finland3)

Is the change of projection wrong? Because when I write Finland3 I get:

class       : SpatialPolygonsDataFrame 
features    : 742 
extent      : 24.29866, 25.20134, 60.23717, 60.68294  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : error, tile_x, tile_y 
min values  :     0,    227,    284 
max values  :     0,    227,    284 

This map shows the tile I want, but the geocoding is wrong at least north-south.

QUESTION:is the projection wrong?????

EDIT AFTER NEW SHP I still get wrong coordinates, this tile has extent (24.29866, 25.20134, 60.23717, 60.68294) (xmin, xmax, ymin, ymax) and if I check these on google maps they are not overlapping, but I need them to overlap my raster when I import I just do:

aa<-readOGR(".", "helsinki_wgs84")
unique(aa@data$tile_x)
unique(aa@data$tile_y)
aa3 <- aa [aa$tile_x == 227 &  aa$tile_y == 284, ]
Finland3<-spTransform(aa3, CRS("+proj=longlat +ellps=WGS84 +no_defs+ellps=WGS84 +towgs84=0,0,0 "))

x11()
plot(Finland3)
  • That is where I got the shapefile. I'll post my code, maybe I am doing something wrong – Irene Jun 30 '14 at 20:25
  • Large Scale Data for Helsinki kartat.kapsi.fi includes coastline datasets – Mapperz Jun 30 '14 at 20:25
  • I am looking at it (it is taking a bit since I don't speak Finnish), in case I upload my code, in case you want to give a look, it is just an IDW interpolation, but I want to crop the raster to the border of the shapefile. – Irene Jun 30 '14 at 20:47
  • Those diva shapes are really poor quality. Can you give the download path to the shape you are using now? I guess it would be much easier to view and cut the data in QGIS. – AndreJ Jul 1 '14 at 16:20
  • @AndreJoost I tried to use those suggested by Kevin R Dyke to no better result (see my answer). I really do not feel confident and miss guidelines, if you could help that would be awesome, really. – Irene Jul 1 '14 at 19:49
3

This set of extracts from OSM data may be what you're looking for. In particular, this shapefile of the coastline around Helsinki.

  • If the extent is to small, take the worldwide version: openstreetmapdata.com/data/coastlines – AndreJ Jul 1 '14 at 10:19
  • The extent is ok, but I have a projection problem, I'll clear up the question. – Irene Jul 1 '14 at 12:17
  • I realised I should have posted another question with this, therefore I did it – Irene Jul 1 '14 at 13:07
  • @Kevin R Dyke can the .prj be possibly wrong? – Irene Jul 1 '14 at 13:49
  • I'm not seeing a problem with the coordinates. Maybe it's the somewhat non-standard projection definition? Try this version, which I quickly made with ogr2ogr. It's the same shapefile, but it's using unprojected WGS84 coordinates. They seem right to me, so it may be worth starting with that and reprojecting as you see fit. – Kevin R Dyke Jul 1 '14 at 15:47
2

At the end I used the shapefile from the National Land survey of Finland even though not perfect, they seemed to be of a higher quality

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.