Почему бы не получить широту и долготу по осям X и Y в функции imshow базовой карты? - PullRequest
Новая
       19

Почему бы не получить широту и долготу по осям X и Y в функции imshow базовой карты?

0 голосов
/ 31 августа 2018

Я использовал широту и долготу в качестве экстента в map.imshow, но не получая широту и долготу по осям x и y, получаю какое-то другое значение v, но в plt.imshow работает нормально, почему ?? Пожалуйста, помогите мне устранить эту проблему благодаря 

from __future__ import print_function

import numpy as np  

import pandas as pd  

import matplotlib.pyplot as plt

from matplotlib_scalebar.scalebar import ScaleBar

from matplotlib_scalebar.scalebar import SI_LENGTH

import matplotlib.font_manager as fm

from mpl_toolkits.basemap import Basemap

import rasterio

from rasterio.plot import show

from matplotlib import colors

from affine import Affine

from pyproj import Proj, transform

#basemap funnction call

map = Basemap(llcrnrlon=-10.5,llcrnrlat=35,urcrnrlon=4.,urcrnrlat=44.,
               projection='tmerc', lat_0 = 39.5, lon_0 = -3.25, 
               suppress_ticks=False)


raster = rasterio.open('lt51350412010304.tif')

#band 3
red = raster.read(3)

red = red.astype(float)

#band 4
nir = raster.read(4)

nir = nir.astype(float)

A = raster.read()

T0 = raster.affine

np.seterr(divide='ignore', invalid='ignore')

#print(raster.shape)
ndvi = np.zeros(raster.shape, dtype=rasterio.float32)

check = np.logical_or ( red > 0, nir > 0 )

ndvi_calc = np.where ( check,  (nir - red ) / ( nir + red ),  0.00000001)

#ndvi calculation 
ndvi = (nir - red ) / ( nir + red )

p1 = Proj(raster.crs)

# All rows and columns
cols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))

# Get affine transform for pixel centres
T1 = T0 * Affine.translation(0.5, 0.5)

# Function to convert pixel row/column index (from 0) to easting/northing at 
centre
rc2en = lambda r, c: (c, r) * T1

Все восточные и северные направления (возможно, есть более быстрый способ сделать это) восточные, северные = np.vectorize (rc2en, otypes = [np.float, np.float]) (строки, столбцы) 

# Project all longitudes, latitudes
p2 = Proj(proj='latlong',datum='WGS84')

longs, lats = transform(p1, p2, eastings, northings)

степень, которую я хочу отобразить 

extents = (longs.min(), longs.max(), lats.min(), lats.max())

print(extents)

#print("crs")
#print(p1)

print('NDVI matrix: ')

print(ndvi_calc)

print('\nMax NDVI: {m}'.format(m=ndvi_calc.max()))

print('Mean NDVI: {m}'.format(m=ndvi_calc.mean()))

print('Median NDVI: {m}'.format(m=np.median(ndvi_calc)))

print('Min NDVI: {m}'.format(m=ndvi_calc.min()))

print('Deviation NDVI: {m}'.format(m=ndvi_calc.std()))

print('Longitude: {m}'.format(m=longs))

print('Latitude: {m}'.format(m=lats))

fig = plt.figure(figsize=(10, 8))

map.drawmapscale(-7., 35.8, -3.25, 39.5, 500, barstyle='fancy')

В той степени, в которой я использую широту и логитус, но получаю что-то другое 

map.imshow(ndvi, aspect='auto', extent = extents)

map.colorbar(location='bottom', pad=0.5, label='Mean amplitude (kA)')

plt.show()
...