'Prisma Satellite Data Projection Problem when exporting as .tiff in Python
I have Prisma Hyperspectral data. However I don't have software to open .he5. That's why I try to convert the data via Python. First of all, I share screenshots about data.
First page of Metadata Second page of Metadata Lat Lon from metadata My WGS84 ENVI UTM ArcGIS UTM to WGS84 Data dimensions
So, My problem is imaging shift (approx. 500m eastwest) on my data (image 4) and ENVI data (image 5 and 7.) I put Image 6 to show that when I convert the data from UTM to WGS84 via ArcGIS, it resample and change dimensions of x and y. Code is:
import h5py
import numpy as np
import rasterio
from osgeo import gdal
input_file='D:\Prisma\PRS_L2D_STD_20200506085209_20200506085213_0001.he5'
f = h5py.File(input_file, 'r')
dataset0 = f.get('HDFEOS')
dataset1 = dataset0.get('SWATHS')
dataset2 = dataset1.get('PRS_L2D_HCO')
dataset3 = dataset2.get('Data Fields')
dataset4 = dataset2.get('Geolocation Fields')
latitude = dataset4.get('Latitude')
longitude= dataset4.get('Longitude')
s_cube = dataset3.get('SWIR_Cube')
v_cube = dataset3.get('VNIR_Cube')
swir = np.array(s_cube)
vnir = np.array(v_cube)
lat = np.array(latitude)
lon = np.array(longitude)
dim=swir.shape
sim=vnir.shape
a=dim[0]
b=dim[2]
c=dim[1]+sim[1]
stacked=np.empty([a,b,c])
i=0
l=0
while i < dim[1]:
stacked[:,:,(238-i)]=swir[:,i,:]
i=i+1
while i < c:
stacked[:,:,(238-i)]=vnir[:,l,:]
i=i+1
l=l+1
stacked = stacked.astype(np.uint16)
stacked_reshaped = np.transpose(stacked, (2, 0, 1))
with rasterio.open( 'deneme.tiff','w',
height=a,
width=b,
count=c,
dtype=stacked_reshaped.dtype,
) as dst:
dst.write(stacked_reshaped)
koordinatsız_tif='C:\\Users\gultekin.erten\Desktop\deneme.tiff'
koordinatlı_tif="C:\\Users\gultekin.erten\Desktop\deneme_son.tiff"
driver = gdal.GetDriverByName('GTiff')
eski_tiff = gdal.Open(koordinatsız_tif)
dtype = gdal.GDT_UInt16
projection = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
bottom_left_y = np.amax(lat)
z=np.amax(lat)-np.amin(lat)
pixel_y_size = z/a
bottom_left_x = np.amin(lon)
p=np.amax(lon)-np.amin(lon)
pixel_x_size = p/b
yeni_tiff = driver.Create(koordinatlı_tif, eski_tiff.RasterXSize,eski_tiff.RasterYSize, 239, dtype)
geotransform = (bottom_left_x, pixel_x_size, 0.0, bottom_left_y,0.0, -pixel_y_size)
yeni_tiff.SetGeoTransform(geotransform)
yeni_tiff.SetProjection(projection)
for j in range(c):
x=eski_tiff.RasterXSize
y=eski_tiff.RasterYSize
array = eski_tiff.GetRasterBand(j+1).ReadAsArray(0,0,x,y)
yeni_tiff.GetRasterBand(j+1).WriteArray(array)
yeni_tiff.FlushCache()
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
Solution | Source |
---|