Revision 8553f079
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 296 | 296 |
lon=raster(paste("HDF4_SDS:UNKNOWN:\"",tswath,"\":1",sep=""))
|
| 297 | 297 |
## HEG Tool reprojection results in large areas of sprious values (regions outside data areas are filled in using the interpolation method) |
| 298 | 298 |
## need to crop the resulting projected data to eliminate these areas |
| 299 |
coords=cbind.data.frame(melt(as.matrix(lat))[,1:3],lon=melt(as.matrix(lon))[,3]) |
|
| 300 |
coords=coords[coords$X1%in%range(coords$X1)|coords$X2%in%range(coords$X2),4:3];colnames(coords)=c("lon","lat")
|
|
| 301 |
#coords=cbind.data.frame(lat=melt(as.matrix(lat))[,3],lon=melt(as.matrix(lon))[,3]) |
|
| 299 |
coords=cbind.data.frame(lat=melt(as.matrix(lat))[,3],lon=melt(as.matrix(lon))[,3],ID=1) |
|
| 302 | 300 |
## crop to big bbox |
| 303 |
# coords=coords[coords$lat<tile_bb$lat_max+0.5&coords$lat>tile_bb$lat_min-0.5&
|
|
| 304 |
# coords$lon>tile_bb$lon_min-0.5&coords$lon<tile_bb$lon_max+0.5,]
|
|
| 305 |
# coordinates(coords)=c("lon","lat")
|
|
| 306 |
# proj4string(coords)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
|
|
| 301 |
coords=coords[coords$lat<tile_bb$lat_max+0.5&coords$lat>tile_bb$lat_min-0.5& |
|
| 302 |
coords$lon>tile_bb$lon_min-0.5&coords$lon<tile_bb$lon_max+0.5,] |
|
| 303 |
coordinates(coords)=c("lon","lat")
|
|
| 304 |
proj4string(coords)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs " |
|
| 307 | 305 |
## project to sinusoidal |
| 308 | 306 |
coords2=spTransform(coords,CRS(projection(td))) |
| 307 |
writeOGR(coords2,dsn=".",layer=sub("[.]hdf","",basename(tswath)),driver="ESRI Shapefile",overwrite=T)
|
|
| 308 |
|
|
| 309 |
system(paste("gdal_grid -ot Byte -a count:radius1=10000:radius2=10000 ",
|
|
| 310 |
" -txe ",paste(bbox(td)[1,],sep="",collapse=" "), |
|
| 311 |
" -tye ",paste(bbox(td)[2,],sep="",collapse=" "), |
|
| 312 |
" -outsize ",paste(td@grid@cells.dim,collapse=" "), |
|
| 313 |
" -l ",sub("[.]hdf","",basename(tswath))," ",
|
|
| 314 |
sub("[.]hdf",".shp",basename(tswath))," mask_",sub("[.]hdf",".tif",basename(tswath)),sep=""))
|
|
| 315 |
|
|
| 316 |
|
|
| 317 |
ps=rasterize(coords2,d2[[2]]) |
|
| 318 |
dist=distance(ps,edge=F) |
|
| 319 |
dist2=dist<10000 |
|
| 320 |
|
|
| 309 | 321 |
## fit alpha hull to draw polygon around region with data |
| 310 |
ah=ahull(coordinates(coords2),alpha=100000) |
|
| 311 |
ah2=ah$x[ah$alpha.extremes,] |
|
| 312 |
ah2=ah$edges[,c("x1","y1")]
|
|
| 322 |
# ah=ahull(coordinates(coords2),alpha=100000)
|
|
| 323 |
# ah2=ah$x[ah$alpha.extremes,]
|
|
| 324 |
# ah2=ah$edges[,c("x1","y1")]
|
|
| 313 | 325 |
|
| 314 |
pp = SpatialPolygons(list(Polygons(list(Polygon(coords[c(1:nrow(coords),1),])),1))) |
|
| 315 |
proj4string(pp)=projection(td) |
|
| 326 |
# pp = SpatialPolygons(list(Polygons(list(Polygon(coords[c(1:nrow(coords),1),])),1)))
|
|
| 327 |
# proj4string(pp)=projection(td)
|
|
| 316 | 328 |
|
| 329 |
|
|
| 330 |
plot(stack(lon,lat)) |
|
| 317 | 331 |
plot(coords,add=F);axis(1);axis(2) |
| 318 |
plot(d2[[2]],add=F)
|
|
| 332 |
plot(d2[[8]],add=F)
|
|
| 319 | 333 |
plot(coords2,add=T);axis(1);axis(2) |
| 320 | 334 |
plot(ah,wpoints=F,add=T,col="red") |
| 321 | 335 |
points(ah2,add=T) |
| ... | ... | |
| 379 | 393 |
raster(readRAST6(paste("SZ_",i,sep="")))
|
| 380 | 394 |
) |
| 381 | 395 |
plot(d2[[2]],add=F) |
| 382 |
plot(coords2,pch=16,cex=.2,add=T)
|
|
| 396 |
points(coords2,pch=16,cex=.2,add=T)
|
|
| 383 | 397 |
plot(pp,add=F,usePolypath = FALSE)#(sp.polygons(pp),usePolypath = FALSE) |
| 398 |
|
|
| 384 | 399 |
} |
| 385 | 400 |
|
| 386 | 401 |
|
Also available in: Unified diff
Adding sensor zenith to quality calculation and improving the filter used to select daily data