标签: raster

看来R中的光栅包不区分GeoTIFF的正旋转和负旋转.我感觉这是因为R忽略了旋转矩阵中的负号.我不够精通,无法深入挖掘raster源代码进行验证,但我确实创建了一个可重现的示例来演示问题:

阅读R徽标并另存为GeoTIFF.

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
proj4string(b) <- crs("+init=epsg:32616")

writeRaster(b, "R.tif")

使用Python向tiff添加旋转

import sys
from osgeo import gdal
from osgeo import osr
import numpy as np
from math import *

def array2TIFF(inputArray,gdalData,datatype,angle,noData,outputTIFF):
#    this script takes a numpy array and saves it to a geotiff
#    given a gdal.Dataset object describing the spatial atributes of the data set
#    the array datatype (as a gdal object) and the name of the output raster, and rotation …

我在C上进行物理实验,Young的干涉实验,我制作了一个打印到file大量像素的程序:

for (i=0; i < width*width; i++)
{
    fwrite(hue(raster_matrix[i]), 1, 3, file);
}

其中hue,当给定值[0..255]时,给出char *3个字节,R,G,B.

我想在我的图像文件中放一个最小的标题,以使这个原始文件成为有效的图像文件.

更简洁:切换自:

offset
0000 : height * width : data } my data, 24bit RGB pixels

至

offset
0000 : dword : magic        \
     : /* ?? */              \
0012 : dword : height         } Header <--> common image file
0016 : dword : width         /
     : /* ?? */             /
0040 : height * …

我想使用类似于下面的R studio(在Arc Map中创建)制作一个情节

我试过以下代码:

# data processing
library(ggplot2)
# spatial
library(raster)
library(rasterVis)
library(rgdal)

#
test <- raster(paste(datafold,'oregon_masked_tmean_2013_12.tif',sep="")) # read the temperature raster
OR<-readOGR(dsn=ORpath, layer="Oregon_10N") # read the Oregon state boundary shapefile

gplot(test) +  
  geom_tile(aes(fill=factor(value),alpha=0.8)) + 
  geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
               fill=NA,color="grey50", size=1)+
  coord_equal()

该代码的输出如下所示:

有几点需要注意.首先,R版本缺少分水岭形状文件.那样就好.

其次,R图中较暗的灰色背景是无数据值.在Arc中,它们不显示,但在R中它们显示为gplot.当我使用光栅包中的"plot"时,它们不显示:

plot(test)

我的问题如下:

1. 如何摆脱'gplot'示例中的深灰色NoData填充？
2. 如何将图例(颜色条)设置为合理(例如在ArcMap和栅格'情节'图例中？)
3. 如何控制色彩图？

要注意,我尝试过很多不同的版本

scale_fill_brewer
scale_fill_manual
scale_fill_gradient

等等,但我得到错误,例如

br <- seq(minValue(test), maxValue(test), len=8)

gplot(test)+
geom_tile(aes(fill=factor(value),alpha=0.8)) +

scale_fill_gradient(breaks = br,labels=sprintf("%.02f", br)) +

geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
             fill=NA,color="grey50", size=1)+
coord_equal()

Regions defined for each Polygons …

我想创建一些GIS图,我想知道R是否可以用于此.以下是一些我在概念上与我想要的类似的情节示例:

• 美国的温度图(或等高线图),颜色(或高度)由州GDP确定.因此,状态边界会在结果图中产生不连续性.
• 美国的温度图,其中海拔高度用于数据.在这种情况下,结果图应在州界之间平滑变化.
• 上述2个图的总和(应用了一些缩放).

我刚刚开始学习R,并想知道它是否适合这种工作.查看ggplot2的coord_map,看起来可能会将数据叠加到美国.但是,获取数据以尊重州界可能非常困难.

有什么建议？

如何从时间序列中提取日期？这是一个时间序列:

x = seq (1, 768)
myts <- ts(x, start=1982, frequency=24)

最初我需要为该rts函数创建一个保存日期/时间数据的向量,观察始于1982年,每月进行2次测量直至2013年.

我正致力于使用Raster和RMarkdown记录一些空间数据的可视化,但我遇到的问题是每个图上方都有一堆负空间.这是RMarkdown代码(稍微简化):

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=12, fig.height=8, echo=FALSE,
                      warning=FALSE, message=FALSE)
```

```{r r-packages}
library(maptools)
library(raster)
library(rgdal)
```
###Description of data
Data are taken from the National Land Cover Database - 2011 and represent land cover at a 30m X 30m resolution.
location of data: [National Land Cover Database - 2011]('http://gisdata.usgs.gov/TDDS/DownloadFile.php?TYPE=nlcd2006&FNAME=nlcd_2006_landcover_2011_edition_2014_10_10.zip')

###Import raster file for US landcover and shapefile for state borders and counties

```{r Import raster file for us landcover}
rfile <- '~/Documents/Data/nlcd_2006_landcover_2011_edition_2014_10_10/nlcd_2006_landcover_2011_edition_2014_10_10.img' #location of raster data
r1 <- raster(rfile)

##Import …

我正在使用R中的"光栅"包从光栅文件创建一些地图.我想创建比较光栅,并排显示几个地图.重要的是,无论每个地图中的值如何,所使用的颜色比例对于所有地图都是相同的.例如,如果地图1的值为0-1,而地图2的值为0-0.5,则值为0.5的单元格在两个地图上应具有相同的颜色.

例如:

• map 1的值从0到1
• map 2的值为0到0.5
• 颜色从红色(最低)到绿色(最高)

我希望值为0.5以在两个地图中具有相同的颜色(即黄色,红色和绿色之间).目前的行为是它在地图1中是黄色,在地图2中是绿色.

我无法找到一种方法来完成这项工作.我看不出有任何方法可以设置用于绘图功能的像素值范围.setMinMax()没有帮助(因为'plot'总是计算值).即使尝试手动设置值(例如g1 @ data @ max < - 10)也不起作用(绘图时会忽略这些值).

最后,制作一堆地图(可能预期在同一色标上绘制所有内容)也不起作用 - 每个地图仍然有自己的色标.

有关如何做到这一点的任何想法？

编辑:

我最终使用的解决方案是:

plot( d, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min(minValue( d )),max(maxValue(d)),length.out=100) ) 

我目前正在尝试将各种"ASC"文件中的信息提取到R中,以便对数据进行分析.

问题是我不确定如何读取文件.我尝试了标准的read.table函数,但所有数字都完全相同(-9999.00).为了排除数据损坏的可能性,我在另一个ASC文件中读取并得到了相同的结果.我唯一知道的是,它们之间的文件大小完全相同.

无论如何我可以阅读这些文件吗？我可以看看任何R包吗？

我绑了这个:

    x = read.table("Dropbox/MVZ/aet2009sep.asc")
    y = read.table("Dropbox/MVZ/aet2009oct.asc")

我的产出是

    > head(x, n =20)
         V1        V2
    1         ncols    3486.0
    2         nrows    4477.0
    3     xllcorner -374495.8
    4     yllcorner -616153.3
    5      cellsize     270.0
    6  NODATA_value   -9999.0
    7      -9999.00   -9999.0
    8      -9999.00   -9999.0
    9      -9999.00   -9999.0
    10     -9999.00   -9999.0
    11     -9999.00   -9999.0
    12     -9999.00   -9999.0
    13     -9999.00   -9999.0
    14     -9999.00   -9999.0
    15     -9999.00   -9999.0
    16     -9999.00   -9999.0
    17     -9999.00   -9999.0
    18     -9999.00   -9999.0
    19     -9999.00   -9999.0
    20     -9999.00   -9999.0

    head(y, n =20) …

我有一组由耦合扫描电子显微镜和能量色散光谱生成的图像和 X 射线数据。这是我的问题：

我像这样对岩石表面的横断面进行了成像（紫色框勾勒出横断面区域）：

我想要真正的高分辨率，所以我使用 7 张放大 3000 倍的图像并将它们与 Photoshop 中的 photomerge 脚本拼接在一起。这是单个图像的示例：

以及它在光合图像横断面上的位置：

在这 7 个位置中的每一个，我还收集了 X 射线数据，为检测到的每个元素生成元素图，并将其写入 TIFF。我还想将每个元素图 TIFF 拼接在一起，以便我可以将其叠加在岩石的合并横断面图像上。这是我想要的结果：

问题是元素贴图没有足够的功能，无法使用 photomerge 将它们拼接在一起。它基本上是二进制的 - 如果检测到元素，则像素是某种颜色（例如，在我的示例图像中，红色代表铁，黄色代表硫），如果未检测到元素，则像素为黑色。您可以看到大部分元素地图大部分是黑色的。

我现在有 ~20 个横断面 x 7 个图像，每个 x ~10 个元素，这导致需要将 ~1400 个图像放在一起，因此需要自动化。

我的想法是用 photomerge 将岩石图像拼接在一起。photomerge 的输出是一个智能对象，其中每个图像都是一个图层。然后，我将使用一个脚本来获取 photomerged 图像对象中 7 个图像中每一个的左上角坐标、宽度和高度。然后，我会将这些属性放置并分配给 7 个图像的每个相应元素映射，以生成“合并”元素映射以覆盖在图像上。我试图自己解决这个问题，但我不精通 javascript 并且无法围绕 Photoshop API。

我在这里在 Github 上上传了一个示例数据集。7 个横断面位置从左到右：-2、-1、0、1、2、3、4。每个位置都有岩石和子目录的图像以及元素数据。

我使用{光栅}夹(或作物)的基础上的不规则shape文件(亚马逊生物群落),但输出总是具有矩形程度的栅格.不过,我需要在shape文件完全相同的几何形状的输出.有小费吗？干杯.

library(raster)
library(rgdal)

myshp <- readOGR("Amazon.shp", layer="Amazon")
e <- extent(myshp)
myraster <- raster("Temperature.tif")
myraster.crop <- crop(myraster, e, snap="out", filename="myoutput.tif")