2

I'm very new to the world of GIS and was working with MOSDAC data from ISRO to extract irradiance parameters such as GHI, DNI and DHI. The output file is in a format called .h5 that I'm not familiar. I tried to open the file in a software called Panoply and found all the data available there. I can also see the data being plotted on a map but for my requirement, I need the final output in lat/long format. Is there anyway to convert the projections in meter to latitude/longitude values as I need this for further processing?

Attaching some screenshots of the data along with the original file.

Link to file: https://drive.google.com/file/d/1xQHNgrlrbyNcb6UyV36xh-7zTfg3f8OQ/view?usp=share_link

enter image description here

enter image description here

2 Answers 2

2

I suppose you want to convert from UTM coordinates to lat and lon. You can do this conveniently with Python. Probably there is a more efficient way/code to do this. But this was the first idea that came to my mind. Remember that I didn't know which EPSG the input data had, so in my example I used an arbitrary EPSG from India, which probably gives wrong values.

import pandas as pd
import numpy as np
import h5py

from osgeo import osr

def convert_coordinates(row, input_epsg):

    output_epsg = 4326
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(input_epsg)
    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(output_epsg)

    coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

    # transform point, height z is optional
    point_x = row['X'] # longitude
    point_y = row['Y'] # latitude
    x, y, z = coordTransform.TransformPoint(point_x, point_y) # could include z variable
    return (x,y)

# read coordinates in a dataframe
df_y = pd.DataFrame(np.array(h5py.File("example.h5")['Y']))
df_x = pd.DataFrame(np.array(h5py.File("example.h5")['X']))

# combine dataframes
df_comb = pd.concat([df_x,df_y], ignore_index=True, sort=False, axis=1)
# name columns
df_comb.set_axis(['X','Y'], axis=1, inplace=True)
# convert coordinates to degree
df_comb[['lat','lon']] = df_comb.apply(convert_coordinates,input_epsg=7755, axis=1, result_type='expand')
# clean up dataframe
df_comb.drop(columns={'X','Y'}, inplace=True)
df_comb.replace([np.inf, -np.inf], np.nan, inplace=True)
df_comb.dropna(inplace=True)

# save result(optional)
df_comb.to_csv("data.csv", index=False)

8
  • I tried running your code with input epsg 3857 as it is mercator and output as 4326. But the lat/long values I'm getting are not matching with the one visualized in panoply. When i tried to visualize the values in google earth, it appears like a straight line across african continent. Nov 30, 2022 at 10:41
  • Yes, strange, I visualized the default X and Y coordinates in Qgis and also got this straight line. Just to clarify: The z-value only describes the intensity of e.g. the DHI, right?
    – zeppeh
    Dec 2, 2022 at 10:58
  • Yeah Z-value has DHI,DHI,DNI,etc based on the dataset. Dec 6, 2022 at 6:51
  • I used this code for conversion: github.com/rishikeshsreehari/boring-stuff-with-python/blob/main/… But still not able to get correct results. Is it a problem with the EPSG codes? Dec 8, 2022 at 10:28
  • Are you sure that the EPSG used is 3857? Otherwise, the codes used should be valid and the conversion correct.
    – zeppeh
    Dec 8, 2022 at 16:52
0

Adding to @zeppeh answer, the following code worked specifically for my requirement. Here the inproj is read from the projection of the file itself.

with h5py.File("mer.h5", "r") as file:
    df_X = pd.DataFrame(file.get("X")[:-2], columns=["X"])
    df_Y = pd.DataFrame(file.get("Y"), columns=["Y"])
    DHI = file.get("DHI")[0][:, :-2].reshape(-1)
    projection= file['Projection_Information']
    lat_0= float(projection.attrs['standard_parallel'])
    lon_0= float(projection.attrs['longitude_of_projection_origin'])
    

final = df_Y.merge(df_X, how="cross").assign(DHI=DHI)[["X", "Y", "DHI"]]


projString='+proj=merc +lat_ts='+str(lat_0)+' +lon_0='+str(lon_0)
inProj=Proj(projString)
outProj=Proj(init='epsg:4326')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.