|
1、准备安装python包apt-getupdateapt-getinstall-ypython3-gdalpipinstallgdalpipinstallpyproj2、下载高程tif数据全国各省12.5m高程数据-免费下载-资源下载-数字地球开放平台注册下,下载他的网盘下载工具,然后选取要使用的区域下载。我下载的宁夏的。3、python代码 以下的代码是专门针对以上的高程数据写的代码。不同的高程数据可能需要不同的投影计算方法。使用Pyproj库进行投影转换: 通过pyproj.Proj方法,确认源投影(EPSG:4326)和目标投影(EPSG:32648)。 将输入的经纬度数据从WGS84转换为UTM48N坐标系统。importgdalimportosrfrompyprojimportProj,transformdefprint_geo_transform_info(transform):print(f"Geotransform:")print(f"Origin{transform[0]},{transform[3]})")print(f"PixelSize{transform[1]},{transform[5]})")defget_elevation(tif_path,lon,lat):#打开GeoTIFF文件dataset=gdal.Open(tif_path)ifdatasetisNone:raiseFileNotFoundError(f"Cannotopenfile:{tif_path}")#打印GeoTIFF文件信息print("GeoTIFFFileInformation:")print(f"Rastersize:{dataset.RasterXSize}x{dataset.RasterYSize}")print(f"Numberofbands:{dataset.RasterCount}")#获取地理变换信息transform_info=dataset.GetGeoTransform()iftransform_infoisNone:raiseValueError("Failedtogetgeotransformfromthedataset.")print_geo_transform_info(transform_info)#获取投影信息proj_info=dataset.GetProjection()ifproj_infoisNone:raiseValueError("Failedtogetprojectionfromthedataset.")print(f"Projection:{proj_info}")#获取栅格波段band=dataset.GetRasterBand(1)ifbandisNone:raiseValueError("Failedtogetrasterbandfromthedataset.")#创建投影转换对象in_proj=Proj(init='epsg:4326')#WGS84out_proj=Proj(init='epsg:32648')#UTMzone48Ntry:#将经纬度坐标转换为UTM坐标srcX,srcY=transform(in_proj,out_proj,lon,lat)print(f"TransformedCoordinates:X={srcX},Y={srcY}")#检查是否在GeoTIFF数据范围内x_min=transform_info[0]x_max=x_min+transform_info[1]*dataset.RasterXSizey_max=transform_info[3]y_min=y_max+transform_info[5]*dataset.RasterYSizeprint(f"GeoTIFFXrange:{x_min}to{x_max}")print(f"GeoTIFFYrange:{y_min}to{y_max}")ifnot(x_min
|
|