Mais

X e Y após a transformação de lat e lon de Albers

X e Y após a transformação de lat e lon de Albers


essa é provavelmente uma pergunta muito óbvia, mas está derretendo meu cérebro:

ao transformar lat e lon em [x, y] (por exemplo, usando este código https://gist.github.com/RandomEtc/476238) Eu obtenho o seguinte: {x: -0,27, y: -1,12}

O que são x e y? Esses são radianos? Como posso usar xey para encontrar o ponto exato para aquele lat e lon no mapa projetado em 2D de maio?

Editar:

Este é o mapa com o qual estou trabalhando: http://www.jetplan.com/weather/data/maps/eutball06.gif ">


Para usar qualquer fórmula, você precisa dos parâmetros da projeção Albers Equal Area:

latitude e longitude de origem primeiro e segundo paralelo

Você pode adivinhar a longitude a partir do único meridiano que é diretamente vertical, mas os outros são mais difíceis de obter.

Você pode se sair melhor usando uma ferramenta de georreferenciamento. O georreferenciador dentro do QGIS faz um bom trabalho.

Para obter o endereço real do pixel, siga as perguntas vinculadas à direita.


Projetar Coordenadas de Latitude-Longitude para x - y

Projete as coordenadas de latitude-longitude para x - y coordenadas especificando uma projeção de mapa. Em seguida, exiba as coordenadas projetadas em um mapa.

Para fazer isso, primeiro especifique as coordenadas de latitude e longitude dos marcos em Boston. Especifique as coordenadas no CRS geográfico NAD83.

Em seguida, importe uma imagem GeoTIFF de Boston como uma matriz e um objeto de referência de células de mapa. Obtenha informações sobre a projeção do mapa consultando a propriedade ProjectedCRS do objeto de referência. Verifique se o CRS geográfico subjacente ao CRS projetado é NAD83.

Projete as coordenadas de latitude-longitude para x - y coordenadas usando o mesmo CRS projetado como a imagem GeoTIFF.

Exibe a imagem GeoTIFF e as coordenadas projetadas no mesmo mapa. Altere o símbolo do marcador e a cor das coordenadas para que fiquem mais visíveis. Em seguida, adicione rótulos de eixo.


GeographicLib

O comentário foi marcado como spam.
Ver e moderar todos os comentários de "Ajuda" postados por este usuário

Eu tenho um sistema que precisa exibir eventos em um mapa, as extensões do mapa são pequenas o suficiente para que eu possa assumir um modelo cartesiano (transformação afim) para os pontos. Mas eu preciso do mapa georreferenciado e meus usuários desejam carregar um GeoTIFF e não fornecer manualmente as coordenadas dos cantos do mapa.

GDAL é capaz de fazer isso diretamente de um arquivo TIFF, ao ler o código-fonte gdalinfo, fui capaz de acessar diretamente as coordenadas do canto do mapa que são impressas por gdalinfo. Mas o GDAL também ocupa uma grande área com muitas dependências e o simples carregamento de sua biblioteca abre as portas do servidor TCP, o que não é aceitável para meu projeto. Portanto, eu prefiro usar GeographicLib para a conversão das informações de georeferência.

Sou capaz de ler automaticamente a estrutura da tag GeoTIFF, criar instâncias das classes de projeção e usar a função Reverse para obter latitude e longitude. Para TransverseMercator, isso está funcionando muito bem, mas para outros esquemas os resultados são significativamente diferentes do que indica GDAL. Eu agradeceria um par de olhos experientes para me ajudar a encontrar a discrepância.

As chamadas que estou fazendo para GeographicLib para este arquivo são as seguintes (C #, valores reais substituídos das variáveis ​​de parâmetro):

a projeção construída é

Como você pode ver, eles variam dos valores de GDAL em cerca de 7 minutos de arco de latitude e 12 minutos de arco de longitude.

Posso estar usando o k errado (escala radial)? É porque o datum [EPSG 6267] não foi contabilizado?


Considerações de compatibilidade

Minvtran será removido

Não recomendado a partir de R2020b

A função minvtran será removida em uma versão futura. Na maioria dos casos, use a função projinv. Se o campo mapprojection dos eixos do mapa atual ou da estrutura de projeção do mapa especificada for 'globo', use a função ecef2geodetic.

Esta tabela mostra alguns usos típicos da função minvtran e como atualizar seu código para usar a função projinv.

Esta tabela mostra alguns usos típicos da função minvtran e como atualizar seu código para usar a função ecef2geodetic.

Se o valor do campo de projeção de mapa da estrutura de projeção de mapa estiver listado nesta tabela, você deve especificar todas as unidades lineares em metros, incluindo coordenadas e campos de estrutura de projeção de mapa. Depois de projetar as coordenadas, você pode convertê-las em outras unidades usando a função unidade de razão.


[x, y] = mfwdtran (lat, lon)
[x, y, z] = mfwdtran (lat, lon, alt)
[. ] = mfwdtran (mstruct.)

[x, y] = mfwdtran (lat, lon) aplica a transformação direta definida pela projeção do mapa nos eixos atuais do mapa. Você pode usar esta função para converter localizações de pontos e vértices de linha e polígono dados em latitudes e longitudes em um sistema de coordenadas de mapa planar projetado.

[x, y, z] = mfwdtran (lat, lon, alt) aplica a projeção direta à entrada 3-D, resultando na saída 3-D. Se a entrada alt estiver vazia ou omitida, então alt = 0 será assumido.

[. ] = mfwdtran (mstruct.) requer uma estrutura de projeção de mapa válida como o primeiro argumento. Neste caso, nenhum eixo de mapa é necessário.


[x, y] = mfwdtran (lat, lon)
[x, y, z] = mfwdtran (lat, lon, alt)
[. ] = mfwdtran (mstruct.)

[x, y] = mfwdtran (lat, lon) aplica a transformação direta definida pela projeção do mapa nos eixos atuais do mapa. Você pode usar esta função para converter localizações de pontos e vértices de linha e polígono dados em latitudes e longitudes em um sistema de coordenadas de mapa planar projetado.

[x, y, z] = mfwdtran (lat, lon, alt) aplica a projeção direta à entrada 3-D, resultando na saída 3-D. Se a entrada alt estiver vazia ou omitida, então alt = 0 será assumido.

[. ] = mfwdtran (mstruct.) requer uma estrutura de projeção de mapa válida como o primeiro argumento. Neste caso, nenhum eixo de mapa é necessário.


Sobreposição de dados geográficos em imagens GeoTIFF

PERGUNTA: Estou tentando exibir uma imagem de dados AVHRR NDVI no formato GeoTIFF em uma projeção de mapa em IDL. Não importa o que eu faça para georregistrar a imagem, os contornos continentais ficam um pouco desfocados. De acordo com a documentação publicada, esta imagem está em uma projeção Albers Conical Equal Area, com os pixels do canto (do canto inferior esquerdo e no sentido horário) localizados nestes locais de latitude / longitude:

Você pode ver na imagem abaixo que os contornos continentais que coloquei em torno do continente africano são quase correto, mas não exatamente correto. O erro é especialmente perceptível em torno de Madegascar, mas você também pode vê-lo ao longo da costa oeste da África. Você pode me ajudar a esclarecer isso?

Esta imagem é quase registrado com a projeção do mapa, mas não exatamente.

RESPOSTA: Como esta é uma imagem GeoTIFF, a própria imagem contém informações que supostamente permitem ao usuário final associar um espaço de coordenadas de dados de projeção de mapa à própria imagem. Os dados, que são chamados metadados, é obtido a partir da imagem GeoTIFF no momento em que a imagem é lida em IDL, fornecendo um GEOTIFF palavra-chave para o READ_TIFF função que lê os dados da imagem. Uma estrutura contendo tags GeoTIFF é retornada pela palavra-chave.

Então você pode jogar junto se quiser, eu disponibilizei uma imagem GeoTIFF para você na minha página web. Extraia a imagem GeoTIFF para o mesmo diretório em que você executará o código IDL de exemplo, que chamei de tiffoverlay.pro.

Para abrir o arquivo GeoTIFF e lê-lo em IDL, digite estes comandos.

Observe que coletamos os metadados GeoTIFF em uma variável de estrutura IDL chamada geotag, e que invertemos a direção Y da imagem. O último é quase sempre necessário ao trabalhar com imagens TIFF em IDL, uma vez que o padrão TIFF define o ponto (0,0) no canto superior esquerdo da imagem, e IDL quase sempre prefere usar o canto inferior esquerdo do imagem para esse fim. Eu sempre inverto os dados reais da imagem, em vez de apenas exibir a imagem "de cabeça para baixo" configurando o !Pedido variável do sistema, visto que muitas vezes desejo interagir com imagens e quero evitar a confusão que é inevitável quando meu local & # 8220 na imagem & # 8221 é diferente de minha & # 8220 exibição da imagem & # 8221 na janela de gráficos.

Os metadados GeoTIFF estão disponíveis como geotags, que são representados em IDL como campos no geotag estrutura.

Infelizmente, a maioria dessas informações provavelmente não significa muito para você, e a documentação do IDL não o ajuda a resolver isso. Para entendê-lo, você provavelmente terá que ir a algum lugar como RemoteSensing.org, onde poderá encontrar documentação sobre o próprio formato GeoTIFF. O software que achei mais útil para entender os arquivos GeoTIFF é Listgeo, que fornece um dump dos metadados GeoTIFF em um formato legível por humanos. Em minhas máquinas Windows, eu uso um invólucro GUI para Listgeo, nomeado listgeo_GUI. Não é um software elegante, mas acho que é essencial.

Para entender as geotags específicas que você encontra em seu arquivo GeoTIFF, você terá que pesquisar as especificações do arquivo GeoTIFF. Normalmente procuro o nome de geotag específico neste documento. As informações de que preciso costumam ser encontradas nos apêndices deste documento.

Aqui está uma parte ligeiramente editada do Listgeo impressão para este arquivo. Isso o ajudará a ver como você pode encontrar as mesmas informações nas tags GeoTIFF diretamente de dentro do IDL.

Vemos que a imagem está em uma projeção de mapa Albers Equal Area. (Podemos obter essas informações dos metadados traduzindo as geotags PROJECTIONGEOKEY e PROJCOORDTRANSGEOKEY - ou campos no geotag estrutura - de acordo com as informações do caderno de encargos do GeoTIFF.) Outras informações essenciais para a configuração da projeção cartográfica também estão disponíveis. Por exemplo, a projeção Albers Equal Area requer dois paralelos padrão e estes são fornecidos nas geotags PROJSTDPARALLEL1GEOKEY e PROJSTDPARALLEL2GEOKEY. O datum usado para esta projeção é WGS84 (uma tradução da geotag GEOGRAPHICTYPEGEOKEY, 4326). Podemos encontrar o centro da projeção do mapa nas geotags PROJNATORIGINLATGEOKEY e PROJNATORIGINLONGGEOKEY.

Além de saber qual projeção de mapa particular foi usada para a imagem, temos que saber onde, exatamente, a imagem está localizada nessa projeção de mapa. Esta informação também está disponível nas geotags, mas você tem que saber analisá-la. De particular importância são os pontos de ligação (MODELTIEPOINTTAG) e a escala da imagem (MODELPIXELSCALETAG). Os pontos de ligação geralmente (mas nem sempre) localizam o canto superior esquerdo da imagem em coordenadas cartesianas. Esta é simplesmente uma superfície 2D plana sobre a qual projetamos a projeção do mapa e sobre a qual desenhamos uma grade em unidades arbitrárias. Nesse caso, as unidades são metros (da tag PROJLINEARUNITSGEOKEY) e a extensão da grade é o diâmetro da Terra. Este sistema de coordenadas arbitrário às vezes é chamado de sistema de coordenadas UV, onde U é a dimensão X (frequentemente associada à longitude) e V é a dimensão Y (freqüentemente associada à latitude).

É neste sistema de coordenadas que os cantos e o centro da imagem são relatados no Listgeo informações acima. Mas podemos calcular os valores de canto de dentro do IDL. Aqui está o código para fazer isso. Neste código, você verá que movi a origem da imagem do canto superior esquerdo para o canto inferior esquerdo. Isso não é estritamente necessário, mas faço isso para minha própria sanidade. Trabalho com IDL há tanto tempo que acabei presumir a origem está no canto inferior esquerdo e eu cometo muitos erros se não estiver. Os pontos finais são encontrados usando os pontos de ligação como coordenadas iniciais e, em seguida, apenas multiplicando o tamanho da imagem em cada dimensão pelo fator de escala para essa dimensão.

Você pode ver na impressão, que este cálculo é idêntico aos cantos relatados por Listgeo.

A próxima etapa é descobrir onde os cantos da imagem estão no espaço de latitude e longitude, ao invés do espaço UV. Para colocar nossa imagem no espaço UV, tivemos que transformar os dados por meio de uma projeção de mapa. (Em outras palavras, nossa imagem representava alguma parte de uma superfície 3D, e nós a transformamos, por meio de uma projeção de mapa, na superfície 2D do sistema de coordenadas UV.) Para voltar a um sistema de coordenadas 3D, representado por latitude e longitude coordenadas, temos que realizar uma transformação inversa. Nós fazemos isso com MAP_PROJ_INVERSE, mas para usá-lo, temos que estabelecer o espaço de projeção do mapa com MAP_PROJ_INIT. O código se parece com isso. (Eu ordenei os argumentos para que as coordenadas comecem no canto inferior esquerdo da imagem e prossigam no sentido horário, para que os valores possam ser comparados com os valores de canto publicados.)

Imprimindo esses valores de latitude / longitude, você pode ver que eles estão não o mesmo que os valores publicados. Depois de passar um bom tempo trabalhando com a documentação de imagens de satélite, não ficaria surpreso se os resultados publicados estivessem simplesmente errados. Isso acontece o tempo todo. Mas, acho que neste caso, os resultados publicados descrevem as localizações do Centro dos pixels do canto, e em IDL precisamos trabalhar com as bordas dos pixels. Portanto, se você quiser usar os valores de canto publicados, terá que mover suas localizações em meio pixel nas dimensões U e V para obter uma sobreposição precisa do mapa.

Para sobrepor dados na projeção do mapa, a extensão da projeção do mapa deve ser exatamente coincidente com a imagem. Assim, temos que conhecer os limites da projeção cartográfica. Podemos tentar calculá-los a partir dos pixels do canto. Escolhendo dois cantos opostos que têm a separação máxima, escrevemos código como este.

Infelizmente, você pode ver (pelos valores calculados de longitude e latitude acima) que o limite não descreve exatamente a localização da imagem no mapa. Realmente precisamos de uma matriz de limite de oito elementos para descrever os quatro cantos da imagem. Algo assim existe em IDL para o MAP_SET comando. (Embora você deva especificar a localização esquerda, superior, direita e inferior da imagem, e não os cantos.) Mas se passarmos uma matriz de oito elementos para o LIMITE palavra-chave em MAP_PROJ_INIT não reclama, mas ignora completamente os valores!

Vamos ver o que acontece se continuarmos a usar o limite conforme calculamos acima. MAP_PROJ_INIT, mas com o LIMITE conjunto de palavras-chave.

O MAP_PROJ_INIT função retorna uma variável de estrutura de mapa. Um dos campos desta variável de estrutura é UV_BOX. Este campo descreve a extensão do mapa em coordenadas UV. Você terá que usar isso UV_BOX para configurar um sistema de coordenadas de dados que você pode desenhar em cima com MAP_CONTINENTS e MAP_GRID.

Primeiro, configuramos a tela e exibimos a imagem. Você precisará de programas da Coyote Library se desejar executar esta parte do código.

Em seguida, configuramos o sistema de coordenadas de dados UV. É aqui que usamos o UV_BOX da estrutura do mapa. Observe o ENREDO O comando é construído de forma que nada seja realmente desenhado na tela. Sua finalidade é simplesmente definir as variáveis ​​de sistema apropriadas que estabelecerão o sistema de coordenadas de dados UV.

Finalmente, podemos traçar limites políticos e linhas de grade. Observe que temos que passar a estrutura do mapa para MAP_CONTINENTS e MAP_GRID portanto, eles têm os meios para converter seus dados de latitude / longitude nas coordenadas UV adequadas.

Você vê o resultado na figura abaixo. Observe que os limites estão ligeiramente deslocados. Você notará isso especialmente em torno de Madagascar, no canto inferior direito da figura.

As sobreposições não bastante combinar. Isso ocorre porque não foi possível especificar a localização da imagem na projeção do mapa com exatidão. Esta é uma limitação da palavra-chave LIMIT de MAP_PROJ_INIT. Ele precisa ser capaz de aceitar uma palavra-chave LIMIT de oito elementos em vez de apenas uma palavra-chave LIMIT de quatro elementos.

Apenas para completar, vamos calcular um elemento de oito LIMITE palavra-chave, como poderíamos fazer se estivéssemos especificando o espaço de coordenadas de dados lat / lon com um MAP_SET comando. O código ficaria assim, no qual estou calculando os lados esquerdo, superior, direito e inferior da imagem.

Você pode ver nos resultados abaixo que o MAP_PROJ_INIT comando leva o elemento de oito LIMITE, mas ignora seus valores.

O comando MAP_PROJ_INIT leva um LIMIT de oito elementos, mas ignora completamente seus valores. O UV_BOX produzido está completamente errado.

Para reiterar, o problema é que não podemos extrair o limite exato da projeção do mapa para esta imagem com um elemento de quatro LIMITE variedade. Mas, na verdade, já temos os limites adequados calculados anteriormente. Estas são as coordenadas que colocamos nas variáveis xOrigin, xEnd, yOrigin, e ano. Simplesmente usamos essas variáveis ​​para configurar nosso espaço de coordenadas de dados. O código se parece com isso.

Você vê os resultados corretos na figura abaixo.

Esta imagem usa os limites calculados do arquivo GeoTiff para configurar o espaço de coordenadas dos dados e o resultado é perfeito.

Uma maneira alternativa de obter resultados perfeitos é usar o UV_BOX a partir de MAP_PROJ_IMAGE via seu UVRANGE palavra-chave, em vez da MAP_PROJ_INIT estrutura do mapa. (Após um ótimo muito esforço, aprendi que é assim imap consegue obter os contornos em torno de seus mapas perfeitos.) MAP_PROJ_IMAGE comando, ao contrário do MAP_PROJ_INIT comando, usa uma palavra-chave de posicionamento de oito elementos para calcular o UV_BOX. Os comandos são semelhantes aos acima. E o resultado parece idêntico à figura diretamente acima.

Os valores de UV_BOX para os quatro métodos são mostrados abaixo. Observe que as coordenadas UV_BOX são idênticas para os dois últimos métodos.

MAP_PROJ_INIT UV_BOX: -4717397,8 -4644644,7 4733971,2 4600353,9 8 elementos UV_BOX: -19170318. -6842903,5 18106165. 7207980,3 UV_BOX calculado: -4599918,3 -4615852,5 4616081,7 4600147,5 MAP_PROJ_IMAGE UV_BOX: -4599918,3 -4615852,5 4616081,7 4600147,5

Muito obrigado a Matt Savoie por sua ajuda na resolução deste mistério de longa data.

Objeto de mapa gráfico Coyote

O programa Coyote Graphics cgGeoMap pode fazer todo esse georregistro para você e criar um objeto de mapa (cgMap) que pode configurar o sistema de coordenadas de dados do mapa para você automaticamente. Ele pode até ler e retornar a imagem para você a partir do arquivo TIFF. Os comandos para fazer isso terão a seguinte aparência.

filename = 'AF03sep15b.n16-VIg.tif' cgDisplay, 500, 500, WID = 5, Title = 'Esboço cgMap Object', $ XPOS = 50, YPOS = 50 alberMap = cgGeoMap (filename, IMAGE = tiffImage, / OnImage) scaled = BytScl (tiffImage, Top = 253) +1 índice = Onde (scaled EQ 1) scaled [index] = 0B TVLCT, cgColor ('ivory', / Triple), 0 cgLoadCT, 33, NColors = 253, Bottom = 1 cgImage, scaled, POSITION = pos, / KEEP_ASPECT Configure um sistema de coordenadas do mapa. alberMap -> Desenhe Desenhe fronteiras e costas políticas. cgMap_Continents, / HIRES, Color = 'Dark Slate Blue', $ Map_Structure = alberMap, / COUNTRIES, / COASTS Adicionar linhas de grade cgMap_Grid, / LABEL, Color = 'Slate Blue', Map_Structure = alberMap

Versão do IDL usada para preparar este artigo: IDL 7.1.

[Retornar para dicas de programação IDL]
Copyright e cópia 2008 David W. Fanning
Escrito em 14 de março de 2008
Última atualização em 22 de novembro de 2011


Exemplos

Antes de usar o minvtran, é necessário criar uma estrutura de projeção do mapa. Você pode fazer isso com axesm ou a função defaultm:

Os seguintes dados de latitude e longitude para o Distrito de Columbia são obtidos do arquivo de forma usastatelo:

Esses dados podem ser projetados em coordenadas cartesianas da projeção de Mercator usando a função projfwd:

Para transformar o projetado x-y dados de volta ao sistema geográfico não projetado, use a função minvtran:


9.1 Preparação de dados em R

  1. Exercício 9.1 - Baixe e extraia a pasta zip em seu local preferido
  2. Defina o diretório de trabalho para a pasta extraída em R em Arquivo - Alterar diretório.
  3. Primeiro, precisamos carregar os pacotes necessários para o exercício

biblioteca (sp)
biblioteca (rede)
biblioteca (rgdal) #readOGR
biblioteca (rgeos) #gIntersection
biblioteca (raster) #para usar a função "raster"
biblioteca (adehabitatHR)
biblioteca (maptools) #readAsciiGrid
biblioteca (zoológico)

nevado & lt-read.csv ("Amostras da neve.csv", cabeçalho = T)
str (nevado)

# Limpe excluindo colunas estranhas, se necessário
nevado e nevado [c (-20: -38, -40: -64)]
nevado $ Status & lt- nevado $ E_schneide

#Faça um quadro de dados espaciais de localizações após remover outliers
coords & lt-data.frame (x = neve $ X_Coordina, y = neve $ Y_Coordina)
utm.crs & lt - "+ proj = utm + zona = 13 + datum = NAD83 + unidades = m + no_defs + ellps = GRS80 + towgs84 = 0,0,0"
utm.spdf & lt- SpatialPointsDataFrame (coords = coords, data = Snow, proj4string = CRS (utm.crs))

#Carregar camada raster DEM
dem & lt-raster ("Snowdem")
imagem (dem)
classe (dem)
proj4string (dem)

#Agora transforme todas as projeções para corresponder ao DEM (ou seja, Albers)
Albers.crs & lt-CRS ("+ proj = aea + lat_1 = 20 + lat_2 = 60 + lat_0 = 40 + lon_0 = -96 + x_0 = 0 + y_0 = 0 + ellps = GRS80 + towgs84 = 0,0,0, 0,0,0,0 + unidades = m + no_defs ")
Snow.spdf & lt-spTransform (utm.spdf, CRS = Albers.crs)

sublette.df & lt- as.data.frame (sublette.spdf)
str (sublette.df)
minx & lt- (min (sublette.df $ x) -10830)
maxx & lt- (max (sublette.df $ x) +10830)
miny & lt- (min (sublette.df $ y) -10830)
maxy & lt- (max (sublette.df $ y) +10830)

## criar vetores dos pontos x e y
x & lt- seq (de = minx, a = maxx, por = 3610)
y & lt- seq (de = miny, a = maxy, por = 3610)

#Código bbox alternativo para pontos espaciais
# mínimo máximo
#x -854784.4 -724665.0
#y 156859.0 247343.2

#Crie vetores dos pontos x e y
#x & lt- seq (de = -854784,4, para = -724665,0, por = 3610)
#y & lt- seq (de = 156859,0, a = 247343,2, por = 3610)

#Crie uma grade de todos os pares de coordenadas (como um data.frame)
xy & lt- expand.grid (x = x, y = y)
classe (xy)
str (xy)

# Identifique a projeção antes de criar um quadro de dados de pontos espaciais
Albers.crs2 & lt - "+ proj = aea + lat_1 = 20 + lat_2 = 60 + lat_0 = 40 + lon_0 = -96 + x_0 = 0 + y_0 = 0 + ellps = GRS80 + towgs84 = 0,0,0,0, 0,0,0 + unidades = m + no_defs "

#NOTA: Albers.crs2 é necessário porque o SPDF precisa de um comando de projeção diferente do spTransform acima

grid.pts & lt-SpatialPointsDataFrame (coords = xy, data = xy, proj4string = CRS (Albers.crs2))
proj4string (grid.pts)
plot (grid.pts)
em grade (grid.pts)
classe (grid.pts)

#Precisa definir pontos como uma grade para converter em polígonos espaciais abaixo
quadriculado (grid.pts) & lt- TRUE
em grade (grid.pts)
str (grid.pts)
plot (grid.pts)

#Converter os pontos da grade em polígonos espaciais, em essência, convertendo-os em um arquivo de forma
gridsp & lt- as (grid.pts, "SpatialPolygons")
str (gridsp)
plot (gridsp)
classe (gridsp)
resumo (gridsp)

grade & lt- SpatialPolygonsDataFrame (gridsp, data = data.frame (id = row.names (gridsp), row.names = row.names (gridsp)))
classe (grade)
plot (grade)
names.grd & lt-sapply (grid @ polygons, function (x) slot (x, "ID"))
texto (coordenadas (grade), rótulos = sapply (slot (grade, "polígonos"), função (i) slot (i, "ID")), cex = 0,3)

# Vamos verificar se todas as células da grade são do mesmo tamanho?
resumo (grade)
getSlots (classe (grade))
classe (slot (grade, "polígonos") [[1]])
getSlots (classe (slot (grade, "polígonos") [[1]]))

#Verifique a área de cada célula na grade de amostragem em metros quadrados
sapply (slot (grade, "polígonos"), função (x) slot (x, "área"))
#[1] 13032100

# Tamanho da célula da grade convertido em quilômetros quadrados
13032100/1000000
# [1] 13.0321 é o tamanho da célula da grade em quilômetros quadrados

# A camada abaixo é uma camada raster HSI de cervo-mula sem perturbação de Sawyer et al. 2009 #incorporado a uma camada porque veados-mula são considerados # hospedeiros do parasita que #estamos investigando

nodis & lt-raster ("snowynodis")
nodis
enredo (nodis)
resumo (nodis)
#Precisa remover NoData da camada HSI de cervo-mula
nodis [is.na (nodis [])] & lt- 0

declive = terreno (dem, opt = 'declive', unidade = 'graus')
aspecto = terreno (dem, opt = 'aspecto', unidade = 'graus')
dem #Agora vamos ver os metadados para cada camada
declive
aspecto

nlcdall & lt- raster ("nlcd_snowy")
nlcdall #Olha os valores raster para a camada de habitat
# Os valores variam de 11 a 95

#Ou traçar para visualizar categorias na legenda
plot (nlcdall)

#Reclassifique os valores em 7 grupos (ou seja, todos os valores entre 0 e 20 são iguais a 1, etc.)
m & lt- c (0, 19, 1, 20, 39, 2, 40, 50, 3, 51,68, 4, 69,79, 5, 80, 88, 6, 89, 99, 7)
rclmat & lt- matrix (m, ncol = 3, byrow = TRUE)
rc & lt- reclassificar (nlcdall, rclmat)
plot (rc) #Agora apenas 7 categorias
rc #Agora apenas 7 categorias
classe (rc)

# Verifique se todos os raster têm a mesma extensão para a criação da pilha
compareRaster (dem, inclinação, aspecto, nodis, rc)
# [1] VERDADEIRO

#############################################################
#############################################################
## NOTA: O código nesta caixa era simplesmente para fins de demonstração para reduzir o tempo geral de processamento durante a aula.
## Ignore esta seção de código se estiver usando seus próprios dados e seu computador tiver os dados adequados
## capacidades de processamento.

#Primeiro, vamos cortar as camadas raster aplicando zoom em apenas alguns locais
plot (rc)
plot (grade, adicionar = T)
pontos (nevado.spdf)

#Code abaixo é usado apenas para aumentar o zoom na grade usando uma camada raster
e & lt- drawExtent ()
#clique no canto superior esquerdo da caixa de corte e no canto inferior direito da caixa de corte para criar zoom
newclip & lt- crop (rc, e)
plot (newclip)
plot (grade, adicionar = T)
pontos (nevado.spdf, col = "vermelho")

#Clip locais dentro da extensão do raster
samp_clip & lt- crop (snowy.spdf, newclip)
plot (newclip)
plot (samp_clip, add = T)
grid_clip & lt- crop (grade, newclip)
plot (grid_clip, add = T)
declive2 & lt- crop (declive, newclip)
aspecto2 & lt- crop (aspecto, newclip)
dem2 & lt- crop (dem, newclip)
HSI & lt- crop (nodis, newclip)

# Verifique se todos os raster têm a mesma extensão para a criação da pilha
compareRaster (dem2, slope2, aspect2, HSI, newclip)
# [1] VERDADEIRO


Sua pergunta realmente não tem nada a ver com “retângulos”, que @MPW já apontou não existem em uma esfera. Existe uma maneira de olhar para o seu problema que o torna um exercício de trigonometria esférica básica, e vou mostrar isso a você. Não afirmo que é a maneira mais rápida de resolver o seu problema.

Você tem um ponto, digamos $ R $, na esfera, digamos que seja Richmond, e você deseja girá-lo sobre um ponto fixo $ O $, digamos que seja Orono. E o seu movimento é uma “rotação” neste sentido: você tem o “título” $ alpha $ de $ O $ a $ R $, que é a bússola-direção que você começa se estiver indo em uma grande -circle path. E você tem a distância $ d $ de $ O $ a $ R $. Sua rotação pergunta a localização que você obterá se percorrer uma distância de $ d $ com um título original de $ alpha + delta $, onde $ delta $ é o ângulo em que você está girando as coisas.

Agora vamos desenhar uma imagem, você terá que fazer o desenho sozinho. Conhecendo a latitude e longitude de Richmond e Orono, você vê que há um triângulo com vértice no pólo norte $ P $. A distância de $ P $ a $ O $ é a colatitude de Orono, que é apenas o complemento da latitude de Orono. Vou chamar isso de $ c_O $. Da mesma forma, você tem a colatitude de Richmond, vou chamar isso de $ c_R $. Então você vê que você tem um triângulo $ RPO $ com as pernas $ c_R $ e $ c_O $, e no pólo, o ângulo é o diferença entre as longitudes das duas cidades. Vou chamar isso de $ lambda $, para manter as letras gregas para os ângulos dos vértices e as letras latinas minúsculas para os comprimentos dos lados, letras maiúsculas para os pontos na esfera.

Então você vê que você tem uma situação SAS, e assim como na trigonometria plana, sua primeira ferramenta a usar é a Lei dos Cossenos, para obter o terceiro lado do seu triângulo, que é a distância de $ O $ a $ R $, que eu ' chamei $ d $. Então, há também uma Lei dos Senos esférica para obter o ângulo do vértice $ alpha $ em $ O = , $ Orono. As leis gerais dos cossenos e senos são $ cos c = cos a cos b + sin a sin b cos gamma ,, $ para cossenos, onde $ gamma $ é o ângulo oposto ao lado do comprimento $ c $. E para a Lei dos Senos, $ frac < sin a> < sin alpha> = frac < sin b> < sin beta> = frac < sin c> < sin gamma> ,, $ onde tenho certeza que você adivinhou que $ beta $ é o ângulo oposto ao lado $ b $ e $ alpha $ é o ângulo oposto ao lado $ a $.

Agora, para nossa configuração de Richmond giratório em torno do centro de Orono: temos nossas duas pernas $ c_R $ e $ c_O $ e nossa diferença de longitude $ lambda $, portanto, obtemos $ cos d = cos c_R cos c_O + sen c_R sin c_O cos lambda ,. $ Agora, com $ d $ e $ lambda $ em mãos, um par de lados e ângulos opostos, você pode usar Sines para obter $ alpha $. Agora adicione seu ângulo de rotação $ delta $ para obter $ alpha '= alpha + delta $, e para encontrar seu Richmond


Assista o vídeo: How to do well in evidence MULTIPLE CHOICE. Be A Better Student