1

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')
prj=ds.GetProjection()
geo_ref=osr.SpatialReference(wkt=prj)

# Reference system for input coordinates
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)

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) 
print(lon,lat)
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.

1 Answer 1

0

The Hydro1k documentation says that the projection is in Lambert Azimuthal Equal Area and the origin of the projection is different for each hydro1k continent file.

hydro1k

Following this post, I used the base LAEA WKT linked and modified the origin coordinates to match the hydro1k documentation for Asia. I was able to get what appears to be the correct projection with the stream.e00 file for asia.

PROJCS["Lambert_Azimuthal_Equal_Area",GEOGCS["Unknown",DATUM["D_unknown",SPHEROID["Unknown",6371007.181,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]]

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.