Skip to content

Instantly share code, notes, and snippets.

@scottyhq

scottyhq/rasterio-vdatum.ipynb

Created April 1, 2022 22:40
Show Gist options
  • Save scottyhq/2153cb01cf951680a4dcb9e1e84b69d7 to your computer and use it in GitHub Desktop.
Save scottyhq/2153cb01cf951680a4dcb9e1e84b69d7 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
rasterio-vdatum.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "ebcc28aa-8e3c-4cdb-94a2-730a4ea57dbe",
"metadata": {},
"source": [
"# Geoid <-> Ellipsoid reprojection\n",
"\n",
"Goal: use rasterio/rioxarray to reproject digital elevation models with Compound CRS or 3D Geodetic CRS"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "380e2c04-ecb8-4d64-8a82-270c711dd8ee",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Rioxarray: 0.10.2\n",
"Rasterio: 1.2.10\n",
"GDAL: 3.4.2\n",
"PROJ: 9.0.0\n"
]
}
],
"source": [
"import rioxarray\n",
"import rasterio\n",
"from osgeo import gdal\n",
"import pyproj\n",
"import logging\n",
"print('Rioxarray:', rioxarray.__version__)\n",
"print('Rasterio:', rasterio.__version__)\n",
"print('GDAL:', gdal.__version__)\n",
"print('PROJ:', pyproj.__proj_version__)"
]
},
{
"cell_type": "markdown",
"id": "70aeaf9d-49ad-45ac-9fbd-075799b337fc",
"metadata": {},
"source": [
"## Copernicus DEM elevations are with respect to GEOID (EPSG:4326+3855) or equivalent (EPSG:9518)\n",
"REF: https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.3\n",
"\n",
"*For elevations with respect to ellipsoid, just need to apply a \"vertical shift grid\":"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "3fde85a1-2e45-4f97-b463-5b583f2a9e53",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Candidate operations found: 1\n",
"-------------------------------------\n",
"Operation No. 1:\n",
"\n",
"unknown id, Inverse of WGS 84 to EGM2008 height (1) + WGS 84 to WGS 84 (G1150), 3 m, World., at least one grid missing\n",
"\n",
"PROJ string:\n",
"+proj=pipeline\n",
" +step +proj=axisswap +order=2,1\n",
" +step +proj=unitconvert +xy_in=deg +xy_out=rad\n",
" +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1\n",
" +step +proj=unitconvert +xy_in=rad +xy_out=deg\n",
" +step +proj=axisswap +order=2,1\n",
"\n",
"Grid us_nga_egm08_25.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/us_nga_egm08_25.tif\n"
]
}
],
"source": [
"!projinfo -s EPSG:4326+3855 -t EPSG:7661 -o PROJ --hide-ballpark --spatial-test intersects"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e6e47757-2d0e-40f9-9aa4-a07f996a8d62",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--2022-04-01 22:39:54-- https://copernicus-dem-30m.s3.eu-central-1.amazonaws.com/Copernicus_DSM_COG_10_N39_00_W109_00_DEM/Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif\n",
"Resolving copernicus-dem-30m.s3.eu-central-1.amazonaws.com (copernicus-dem-30m.s3.eu-central-1.amazonaws.com)... 52.219.74.9\n",
"Connecting to copernicus-dem-30m.s3.eu-central-1.amazonaws.com (copernicus-dem-30m.s3.eu-central-1.amazonaws.com)|52.219.74.9|:443... connected.\n",
"HTTP request sent, awaiting response... 200 OK\n",
"Length: 40712353 (39M) [image/tiff]\n",
"Saving to: ‘Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif’\n",
"\n",
"Copernicus_DSM_COG_ 100%[===================>] 38.83M 10.2MB/s in 3.8s \n",
"\n",
"2022-04-01 22:39:58 (10.2 MB/s) - ‘Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif’ saved [40712353/40712353]\n",
"\n"
]
}
],
"source": [
"# Grab a local copy\n",
"!wget https://copernicus-dem-30m.s3.eu-central-1.amazonaws.com/Copernicus_DSM_COG_10_N39_00_W109_00_DEM/Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "0d1df1c2-eb0b-4901-8eef-cf8c4197870c",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"PROJ: pj_open_lib(proj.db): call fopen(/srv/conda/envs/notebook/share/proj/proj.db) - succeeded\n",
"GDAL: GDALOpen(Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif, this=0x565017c63140) succeeds as GTiff.\n",
"GDAL: Using GTiff driver\n",
"PROJ: pj_open_lib(proj.ini): call fopen(/srv/conda/envs/notebook/share/proj/proj.ini) - succeeded\n",
"GDAL: Computing area of interest: -109, 39.0001, -108, 40.0001\n",
"OGRCT: Wrap source at -108.5.\n",
"PROJ: pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed\n",
"PROJ: pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed\n",
"PROJ: Using https://cdn.proj.org/us_nga_egm08_25.tif\n",
"PROJ: pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed\n",
"PROJ: pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed\n",
"PROJ: Using https://cdn.proj.org/us_nga_egm08_25.tif\n",
"GDAL: GDALDriver::Create(GTiff,./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid.tif,3600,3600,1,Float32,(nil))\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Creating output file that is 3600P x 3600L.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARP: Copying metadata from first source to destination dataset\n",
"GTiff: ScanDirectories()\n",
"GTiff: Opened 1800x1800 overview.\n",
"GTiff: Opened 900x900 overview.\n",
"GTiff: Opened 450x450 overview.\n",
"GDALWARP: Defining SKIP_NOSOURCE=YES\n",
"PROJ: Using coordinate operation WGS 84 to EGM2008 height (1)\n",
"PROJ: pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed\n",
"PROJ: pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed\n",
"PROJ: Using https://cdn.proj.org/us_nga_egm08_25.tif\n",
"PROJ: Using coordinate operation WGS 84 to EGM2008 height (1)\n",
"PROJ: pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed\n",
"PROJ: pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed\n",
"PROJ: Using https://cdn.proj.org/us_nga_egm08_25.tif\n",
"GDAL: GDAL_CACHEMAX = 387 MB\n",
"GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyFloat() Src=0,0,3600x1800 Dst=0,0,3600x1800\n",
"GDAL: GDALWarpKernel()::GWKNearestNoMasksOrDstDensityOnlyFloat() Src=0,1800,3600x1800 Dst=0,1800,3600x1800\n",
"GDAL: Flushing dirty blocks: 0...10...20...30...40...50...60...70...80...90...100 - done.\n",
"GDAL: GDALClose(./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid.tif, this=0x565017fe6c10)\n",
"GDAL: GDALClose(Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif, this=0x565017c63140)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Processing Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.\n"
]
}
],
"source": [
"%%bash\n",
"\n",
"# GDAL knows how to go fetch a vshift grid from the web and apply it:\n",
"\n",
"#INPUT='/vsis3/copernicus-dem-30m/Copernicus_DSM_COG_10_N39_00_W109_00_DEM/Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif' \n",
"INPUT='Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif'\n",
"OUTPUT='./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid.tif' \n",
"\n",
"AWS_NO_SIGN_REQUEST=YES \\\n",
"GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \\\n",
"CPL_DEBUG=ON \\\n",
"PROJ_DEBUG=2 \\\n",
"gdalwarp -s_srs EPSG:4326+3855 -t_srs EPSG:7661 $INPUT $OUTPUT"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "450ac2b2-b689-4666-96ef-556904e3228f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"< STATISTICS_MAXIMUM=3327.7126464844\n",
"< STATISTICS_MEAN=2033.3877791249\n",
"< STATISTICS_MINIMUM=1324\n",
"< STATISTICS_STDDEV=372.76572182465\n",
"> STATISTICS_MAXIMUM=3312.0522460938\n",
"> STATISTICS_MEAN=2016.9418973372\n",
"> STATISTICS_MINIMUM=1306.1044921875\n",
"> STATISTICS_STDDEV=373.10998217841\n"
]
}
],
"source": [
"%%bash \n",
"\n",
"gdalinfo -stats Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif > geoid.info\n",
"gdalinfo -stats Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid.tif > ellipsoid.info\n",
"diff geoid.info ellipsoid.info | grep STATISTICS"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d7f2a5b4-69b4-42fe-868e-1b1fd0d30c78",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"DEBUG:rasterio.env:Entering env context: <rasterio.env.Env object at 0x7f6cb6fbdf10>\n",
"DEBUG:rasterio.env:Starting outermost env\n",
"DEBUG:rasterio.env:No GDAL environment exists\n",
"DEBUG:rasterio.env:New GDAL environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> created\n",
"DEBUG:rasterio._env:GDAL_DATA found in environment.\n",
"DEBUG:rasterio._env:PROJ_LIB found in environment.\n",
"DEBUG:rasterio._env:Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:Entered env context: <rasterio.env.Env object at 0x7f6cb6fbdf10>\n",
"DEBUG:rasterio.env:Got a copy of environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio.env:Entering env context: <rasterio.env.Env object at 0x7f6cb6fbbbe0>\n",
"DEBUG:rasterio.env:Got a copy of environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio.env:Entered env context: <rasterio.env.Env object at 0x7f6cb6fbbbe0>\n",
"DEBUG:rasterio._base:Sharing flag: 0\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALOpen(./Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif, this=0x5633d6e830a0) succeeds as GTiff.\n",
"DEBUG:rasterio._base:Nodata success: 0, Nodata value: -10000000000.000000\n",
"DEBUG:rasterio._base:Dataset <open DatasetReader name='./Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif' mode='r'> is started.\n",
"DEBUG:rasterio.env:Exiting env context: <rasterio.env.Env object at 0x7f6cb6fbbbe0>\n",
"DEBUG:rasterio.env:Cleared existing <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio._env:Stopped GDALEnv <rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:No GDAL environment exists\n",
"DEBUG:rasterio.env:New GDAL environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> created\n",
"DEBUG:rasterio._env:GDAL_DATA found in environment.\n",
"DEBUG:rasterio._env:PROJ_LIB found in environment.\n",
"DEBUG:rasterio._env:Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:Exited env context: <rasterio.env.Env object at 0x7f6cb6fbbbe0>\n",
"DEBUG:rasterio._env:CPLE_None in GTiff: ScanDirectories()\n",
"DEBUG:rasterio._env:CPLE_None in GTiff: Opened 1800x1800 overview.\n",
"DEBUG:rasterio._env:CPLE_None in GTiff: Opened 900x900 overview.\n",
"DEBUG:rasterio._env:CPLE_None in GTiff: Opened 450x450 overview.\n",
"DEBUG:rasterio._crs:Matched. confidence=100, c_code=b'4326', c_name=b'EPSG'\n",
"DEBUG:rasterio._crs:Matched. confidence=100, c_code=b'9518', c_name=b'EPSG'\n",
"DEBUG:rasterio._env:CPLE_None in VRT: No valid sources found for band in VRT file \n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALOpen(<VRTDataset rasterYSize=\"3600\" rasterXSize=\"3600\"><VRTRasterBand /><SRS>COMPD_CS[\"WGS 84 + EGM2008 height\",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.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]],VERT_CS[\"EGM2008 height\",VERT_DATUM[\"EGM2008 geoid\",2005,AUTHORITY[\"EPSG\",\"1027\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"3855\"]],AUTHORITY[\"EPSG\",\"9518\"]]</SRS><GeoTransform>-109.00013888888888,0.0002777777777777738,0.0,40.00013888888889,0.0,-0.00027777777777777973</GeoTransform></VRTDataset>, this=0x5633d6e0aad0) succeeds as VRT.\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: Computing area of interest: -109, 39.0001, -108, 40.0001\n",
"DEBUG:rasterio._env:CPLE_None in OGRCT: Wrap source at -108.5.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"EPSG:4326\n",
"EPSG:9518\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"DEBUG:rasterio._warp:Created transformer and warp output.\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALClose(<VRTDataset rasterYSize=\"3600\" rasterXSize=\"3600\"><VRTRasterBand /><SRS>COMPD_CS[\"WGS 84 + EGM2008 height\",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.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]],VERT_CS[\"EGM2008 height\",VERT_DATUM[\"EGM2008 geoid\",2005,AUTHORITY[\"EPSG\",\"1027\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"3855\"]],AUTHORITY[\"EPSG\",\"9518\"]]</SRS><GeoTransform>-109.00013888888888,0.0002777777777777738,0.0,40.00013888888889,0.0,-0.00027777777777777973</GeoTransform></VRTDataset>, this=0x5633d6e0aad0)\n",
"DEBUG:rasterio._io:Output nodata value read from file: None\n",
"DEBUG:rasterio._io:Output nodata values: [None]\n",
"DEBUG:rasterio._io:all_valid: True\n",
"DEBUG:rasterio._io:mask_flags: ([<MaskFlags.all_valid: 1>],)\n",
"DEBUG:rasterio._io:Jump straight to _read()\n",
"DEBUG:rasterio._io:Window: Window(col_off=0, row_off=0, width=3600, height=3600)\n",
"DEBUG:rasterio._io:IO window xoff=0.0 yoff=0.0 width=3600.0 height=3600.0\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDAL_CACHEMAX = 387 MB\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALDriver::Create(MEM,771c1832-6816-4805-bcaa-2227d049ce6d,3600,3600,1,Float32,(nil))\n",
"DEBUG:rasterio._crs:Matched. confidence=100, c_code=b'9518', c_name=b'EPSG'\n",
"DEBUG:rasterio._io:Set CRS on temp dataset: b'COMPD_CS[\"WGS 84 + EGM2008 height\",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.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]],VERT_CS[\"EGM2008 height\",VERT_DATUM[\"EGM2008 geoid\",2005,AUTHORITY[\"EPSG\",\"1027\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"3855\"]],AUTHORITY[\"EPSG\",\"9518\"]]'\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALDriver::Create(MEM,6930c86a-e1a4-4304-b6cc-fe2cff610e70,3600,3600,1,Float32,(nil))\n",
"DEBUG:rasterio._io:Set CRS on temp dataset: b'GEOGCRS[\"WGS 84 (G1150)\",DYNAMIC[FRAMEEPOCH[2001]],DATUM[\"World Geodetic System 1984 (G1150)\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[ellipsoidal,3],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"ellipsoidal height (h)\",up,ORDER[3],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Geodesy. Navigation and positioning using GPS satellite system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",7661]]'\n",
"DEBUG:rasterio._warp:Created temp destination dataset.\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: Computing area of interest: -109, 39.0001, -108, 40.0001\n",
"DEBUG:rasterio._warp:Created approximate transformer\n",
"DEBUG:rasterio._warp:Created transformer and options.\n",
"DEBUG:rasterio._warp:Setting NUM_THREADS option: 1\n",
"DEBUG:rasterio._warp:Configured to warp src band 1 to destination band 1\n",
"DEBUG:rasterio._warp:Set transformer options\n",
"DEBUG:rasterio._warp:Chunk and warp window: 0, 0, 3600, 3600.\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALWarpKernel()::GWKRealCase() Src=0,0,1800x1800 Dst=0,0,1800x1800\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALWarpKernel()::GWKRealCase() Src=1800,0,1800x1800 Dst=1800,0,1800x1800\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALWarpKernel()::GWKRealCase() Src=0,1800,1800x1800 Dst=0,1800,1800x1800\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALWarpKernel()::GWKRealCase() Src=1800,1800,1800x1800 Dst=1800,1800,1800x1800\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALClose(Temporary destination dataset for _reproject(), this=0x5633d70cb510)\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALClose(771c1832-6816-4805-bcaa-2227d049ce6d, this=0x5633d5072250)\n",
"DEBUG:rasterio._crs:Matched. confidence=100, c_code=b'7661', c_name=b'EPSG'\n",
"/srv/conda/envs/notebook/lib/python3.9/site-packages/rioxarray/raster_writer.py:108: UserWarning: The nodata value (3.402823466e+38) has been automatically changed to (3.4028234663852886e+38) to match the dtype of the data.\n",
" warnings.warn(\n",
"DEBUG:rasterio.env:Got a copy of environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio.env:Entering env context: <rasterio.env.Env object at 0x7f6cb65a0f70>\n",
"DEBUG:rasterio.env:Got a copy of environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio.env:Entered env context: <rasterio.env.Env object at 0x7f6cb65a0f70>\n",
"DEBUG:rasterio._io:Path: UnparsedPath(path='./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif'), mode: w, driver: GTiff\n",
"DEBUG:rasterio._io:Skipped delete for overwrite. Dataset does not exist: './Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif'\n",
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALDriver::Create(GTiff,./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif,3600,3600,1,Float32,(nil))\n",
"DEBUG:rasterio._base:Nodata success: 1, Nodata value: 340282346638528859811704183484516925440.000000\n",
"DEBUG:rasterio.env:Exiting env context: <rasterio.env.Env object at 0x7f6cb65a0f70>\n",
"DEBUG:rasterio.env:Cleared existing <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio._env:Stopped GDALEnv <rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:No GDAL environment exists\n",
"DEBUG:rasterio.env:New GDAL environment <rasterio._env.GDALEnv object at 0x7f6cb785d790> created\n",
"DEBUG:rasterio._env:GDAL_DATA found in environment.\n",
"DEBUG:rasterio._env:PROJ_LIB found in environment.\n",
"DEBUG:rasterio._env:Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:Exited env context: <rasterio.env.Env object at 0x7f6cb65a0f70>\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"EPSG:7661\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"DEBUG:rasterio._env:CPLE_None in GDAL: GDALClose(./Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif, this=0x5633d72b7960)\n",
"DEBUG:rasterio.env:Exiting env context: <rasterio.env.Env object at 0x7f6cb6fbdf10>\n",
"DEBUG:rasterio.env:Cleared existing <rasterio._env.GDALEnv object at 0x7f6cb785d790> options\n",
"DEBUG:rasterio._env:Stopped GDALEnv <rasterio._env.GDALEnv object at 0x7f6cb785d790>.\n",
"DEBUG:rasterio.env:Exiting outermost env\n",
"DEBUG:rasterio.env:Exited env context: <rasterio.env.Env object at 0x7f6cb6fbdf10>\n"
]
}
],
"source": [
"# Rioxarray / rasterio don't seem to handle 3D CRS\n",
"\n",
"logging.basicConfig(level=logging.DEBUG,\n",
" #format='%(asctime)s %(message)s',\n",
" handlers=[logging.StreamHandler(),\n",
" logging.FileHandler('rioxarray.log')], #NOTE: appends by default\n",
" )\n",
"\n",
"\n",
"INPUT = './Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif'\n",
"OUTPUT= './Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif'\n",
"\n",
"with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',\n",
" AWS_NO_SIGN_REQUEST='YES',\n",
" CPL_DEBUG=True,\n",
" PROJ_DEBUG=2, #NOT sure how to expose these\n",
" ):\n",
"\n",
" da = rioxarray.open_rasterio(INPUT)\n",
" print(da.rio.crs)\n",
" \n",
" # Assign true 3D geodetic CRS\n",
" CRS = da.rio.crs.from_epsg(9518) #shorthand for \"EPSG:4326+3855\"\n",
" daGeoid = da.rio.write_crs(CRS)\n",
" print(daGeoid.rio.crs)\n",
" \n",
" # Reproject to Ellipsoid (apply vshift grid)\n",
" daEllipsoid = daGeoid.rio.reproject(7661)\n",
" print(daEllipsoid.rio.crs)\n",
" \n",
" daEllipsoid.rio.to_raster(OUTPUT)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "56687342-6a70-48b7-874b-eba959369071",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"< STATISTICS_MEAN=2033.3877791249\n",
"> STATISTICS_MEAN=2033.3877791246\n",
"< STATISTICS_STDDEV=372.76572182465\n",
"> STATISTICS_STDDEV=372.76572182466\n"
]
}
],
"source": [
"%%bash \n",
"\n",
"gdalinfo -stats Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif > geoid.info\n",
"gdalinfo -stats Copernicus_DSM_COG_10_N39_00_W109_00_DEM_ellipsoid_rio.tif > ellipsoid.info\n",
"diff geoid.info ellipsoid.info | grep STATISTICS"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment