The Hydro1k database for Asia is provided by USGS. It is available here as GTOPO30 HYDRO 1K.

I tried to read the slope values for the Kumamoto area in Japan but the results all came out as 55537(same as for the sea). The slope is shown different when I open it in QGIS.

Here is the prj file:

PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["D_Sphere_Clarke_1866_Authalic",SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

Here is my code:

import osgeo
from osgeo import gdal,osr
import ogr
geo = gdal.Open('/home/Downloads/gt30h1kas/as_slope.bil') #slope file

# Read the projection from prj file
ds = gdal.Open('/home/Downloads/gt30h1kas/as_dem.prj')

# Reference system for input coordinates

ctran2=ogr.osr.CoordinateTransformation(point_ref, geo_ref)

arr = geo.ReadAsArray()

geotransform = geo.GetGeoTransform() 
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

# Kumamoto (32.8031° N, 130.7079° E)
[lon,lat,z]=ctran2.TransformPoint(130.7079, 32.8031) 
offx = int( (lon-originX)/pixelWidth )
offy = int( (lat-originY)/pixelHeight )
print(offx, offy, arr[offx][offy])

Here is the output:

2810038.798891878 -837540.6774932265
7165 7281 55537

Please help me find my mistake.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.