Skip to content

Instantly share code, notes, and snippets.

@thareUSGS

thareUSGS/02_projection_tips.ipynb

Last active June 24, 2020 23:04
Show Gist options
  • Save thareUSGS/8f29b897847a58ba2a0060cfd17adf83 to your computer and use it in GitHub Desktop.
Save thareUSGS/8f29b897847a58ba2a0060cfd17adf83 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
02_projection_tips.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# GDAL Map Projection Tips"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Requirements:\n",
"> conda install -c conda-forge gdal=3 jupyter matplotlib tuiview rasterio\n",
"\n",
"**Data for lesson**\n",
"\n",
"TES Image: \n",
"\n",
"https://astrogeology.usgs.gov/search/map/Mars/GlobalSurveyor/TES/Mars_MGS_TES_Albedo_mosaic_global_7410m\n",
"\n",
"This file is in GeoTIFF format with optional ISIS3 and PDS3 labels available for download. Copy into a \"tes/\" directory below this notebook.\n",
"\n",
"files:\n",
"\n",
"* Mars_MGS_TES_Albedo_mosaic_global_7410m.tif\n",
"\n",
"---\n",
"\n",
"Note: \n",
"* Using **!** allows one to run a routine outside of the notebook in the terminal\n",
"* $variable, sends the variable's value to the terminal (command-line)\n",
"* \"//\" for path works in Linux (Mac) and Windows"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gathering information from GDAL"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here the overall goal is to learn how to crop using a degree boundary so that we can map project to a south Polar Sterographic map projection. The lesson here might seem long-winded, but we are really just learing different ways to grab information from images.\n",
"\n",
"GDAL really likes to work in the current units of the map projection which is generally meters. This can be hard to understand how to crop using latitude and logitudes bounds.\n",
"\n",
"For example, using command-line gdalinfo you can see the horizontal extents in meters for this sample Mars Albedo image Mars_MGS_TES_Albedo_mosaic_global_7410m.tif"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"GDAL's version is: 3.0.4\n"
]
}
],
"source": [
"# First, from GDAL's Python API, let's check GDAL's version.\n",
"# This will also double check our GDAL Python environment is working.\n",
"\n",
"# Import the \"gdal\" submodule from within the \"osgeo\" module\n",
"from osgeo import gdal\n",
"\n",
"# We can check which version we're running by printing the \"__version__\" variable\n",
"print(\"GDAL's version is: \" + gdal.__version__)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# define test image and output file\n",
"global_img = \"tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif\"\n",
"out_SPole_img = \"tes//out_south_pole.tif\""
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# check out the test image\n",
"!tuiview $global_img"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[' \"cornerCoordinates\":{',\n",
" ' \"upperLeft\":[',\n",
" ' -10670400.0,',\n",
" ' 5335200.0',\n",
" ' ],',\n",
" ' \"lowerLeft\":[',\n",
" ' -10670400.0,',\n",
" ' -5335200.0',\n",
" ' ],',\n",
" ' \"lowerRight\":[',\n",
" ' 10670400.0,',\n",
" ' -5335200.0',\n",
" ' ],',\n",
" ' \"upperRight\":[',\n",
" ' 10670400.0,',\n",
" ' 5335200.0',\n",
" ' ],',\n",
" ' \"center\":[',\n",
" ' 0.0,',\n",
" ' 0.0',\n",
" ' ]']\n"
]
}
],
"source": [
"# Just to introduce gdalinfo from the command-line \n",
"# using the \"-json\" flag and pulling out a section\n",
"import pprint\n",
"\n",
"img_info = !gdalinfo -json $global_img\n",
"index = img_info.index(' \"cornerCoordinates\":{')\n",
"#pprint.pprint(img_info)\n",
"pprint.pprint(img_info[index:index + 21])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So using the command-line works, put for scripting we should really use the GDAL **Python API**\n",
"from: https://ceholden.github.io/open-geo-tutorial/python/chapter_1_GDALDataset.html"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{'center': [0.0, 0.0],\n",
" 'lowerLeft': [-10670400.0, -5335200.0],\n",
" 'lowerRight': [10670400.0, -5335200.0],\n",
" 'upperLeft': [-10670400.0, 5335200.0],\n",
" 'upperRight': [10670400.0, 5335200.0]}\n",
"\n",
"{'coordinates': [[[179.9738263, 90.0130868],\n",
" [179.9738263, -90.0130868],\n",
" [-179.9738263, -90.0130868],\n",
" [-179.9738263, 90.0130868],\n",
" [179.9738263, 90.0130868]]],\n",
" 'type': 'Polygon'}\n"
]
}
],
"source": [
"import json\n",
"import pprint\n",
"\n",
"# open a original image dataset\n",
"dataset = gdal.Open(global_img, gdal.GA_ReadOnly)\n",
"json_data = gdal.Info(dataset, format = 'json')\n",
"#type(json_data)\n",
"\n",
"json_formatted_str = json.dumps(json_data, indent=2)\n",
"#pprint.pprint(json_formatted_str)\n",
"\n",
"#close the open handle\n",
"dataset = None\n",
"\n",
"# grab sections from the JSON dictionary and print to notebook\n",
"cornerCoordinates = json_data[\"cornerCoordinates\"]\n",
"pprint.pprint(cornerCoordinates)\n",
"\n",
"print()\n",
"\n",
"coordinates = json_data[\"extent\"]\n",
"pprint.pprint(coordinates)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But this still requires awkward parsing to get individual values, so now just use built-in methods\n",
"\n",
"from: https://download.osgeo.org/gdal/workshop/foss4ge2015/workshop_gdal.html\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Width: 2880\n",
"Height: 1440\n",
"Number of bands: 1\n",
"Pixel size: (7410.0, -7410.0)\n",
"Center: [0.0, 0.0]\n",
"Origin: (-10670400.0, 5335200.0)\n",
"Upper Left Corner: [-10670400.0, 5335200.0]\n",
"Upper Right Corner: [10670400.0, 5335200.0]\n",
"Lower Left Corner: [-10670400.0, -5335200.0]\n",
"Lower Right Corner: [10670400.0, -5335200.0]\n"
]
}
],
"source": [
"from osgeo import gdal\n",
"ds = gdal.Open(global_img)\n",
"print('Width:', ds.RasterXSize)\n",
"print('Height:', ds.RasterYSize)\n",
"print('Number of bands:', ds.RasterCount)\n",
"gt = ds.GetGeoTransform() # captures origin and pixel size\n",
"print('Pixel size:', (gt[1], gt[5]))\n",
"print('Center:', gdal.ApplyGeoTransform(gt,ds.RasterXSize/2,ds.RasterYSize/2))\n",
"print('Origin:', (gt[0], gt[3]))\n",
"print('Upper Left Corner:', gdal.ApplyGeoTransform(gt,0,0))\n",
"print('Upper Right Corner:', gdal.ApplyGeoTransform(gt,ds.RasterXSize,0))\n",
"print('Lower Left Corner:', gdal.ApplyGeoTransform(gt,0,ds.RasterYSize))\n",
"print('Lower Right Corner:',gdal.ApplyGeoTransform(gt,ds.RasterXSize,ds.RasterYSize))\n",
"ds = None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### how to show coordinates in different forms (meters and degrees)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are a couple GDAL utilities to convert from Lon/Lat to meters (or vice-versa). One specific program is *gdaltransform*. To use *gdaltransform* you need an input and out projection. To get the to and from projections,follow steps (1) and (2).\n",
"\n",
"(1) The input projection can be found using the input file.\n",
"\n",
"(2) For the output \"projection\", we want to report lat/lon thus we a string which simply defines Mars. To see a Mars definition click on the Proj4 link from: https://spatialreference.org/ref/iau2000/49900/"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs\n"
]
}
],
"source": [
"# (1.a) input projection using command-line for the input image\n",
"# This method is not recommended, but shows how to simply grab a projection\n",
"# string from the command-line (which can be used directly in gdaltransform)\n",
"srs = !gdalsrsinfo -o proj4 $global_img\n",
"srs = srs[1]\n",
"print(srs)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Image projection (Well Known Text):\n",
"\n",
"PROJCS[\"SimpleCylindrical Mars\",GEOGCS[\"GCS_Mars\",DATUM[\"D_Mars\",SPHEROID[\"Mars\",3396000,0]],PRIMEM[\"Reference_Meridian\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Equirectangular\"],PARAMETER[\"standard_parallel_1\",0],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]\n",
"\n",
"Image projection (Proj4):\n",
"\n",
"+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs\n"
]
}
],
"source": [
"# (1.b) now using Python\n",
"\n",
"# open a original image dataset\n",
"dataset = gdal.Open(global_img, gdal.GA_ReadOnly)\n",
"wkt_proj = dataset.GetProjection()\n",
"print('Image projection (Well Known Text):\\n')\n",
"print(wkt_proj + '\\n')\n",
"\n",
"from osgeo import osr\n",
"inSRS = osr.SpatialReference() # makes an empty spatial ref object\n",
"inSRS.ImportFromWkt(wkt_proj) # populates the spatial ref object with our WKT SRS\n",
"proj4 = inSRS.ExportToProj4() # Exports an SRS ref as a Proj4 string usable by PyProj\n",
"print('Image projection (Proj4):\\n')\n",
"print(proj4)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PROJCS[\"SimpleCylindrical Mars\",GEOGCS[\"GCS_Mars\",DATUM[\"D_Mars\",SPHEROID[\"Mars\",3396000,0]],PRIMEM[\"Reference_Meridian\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Equirectangular\"],PARAMETER[\"standard_parallel_1\",0],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]\n",
"\n",
"+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396000 +units=m +no_defs=True\n"
]
}
],
"source": [
"# Bonus - now use RasterIO Python, nice and easy.\n",
"\n",
"import rasterio\n",
"src = rasterio.open(global_img)\n",
"print(src.crs.wkt + \"\\n\")\n",
"\n",
"srcp4 = src.crs.to_proj4()\n",
"print(srcp4)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# (2) Lets make a proj4 string for a spherical Mars (in degrees).\n",
"\n",
"# Note the radius used in the TES Albedo map in not the recommended Mars radius\n",
"# as defined by the IAU. We will use the IAU recommended radius later.\n",
"\n",
"Mars_truncatedSphere = \"+proj=longlat +a=3396000 +b=3396000 +no_defs\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can find out what Lon/Lat maps into X/Y meters in the image's map projection"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Usage: gdaltransform [--help-general]\n",
" [-i] [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]\n",
" [-ct proj_string] [-order n] [-tps] [-rpc] [-geoloc] \n",
" [-gcp pixel line easting northing [elevation]]* [-output_xy]\n",
" [srcfile [dstfile]]\n",
"\n",
"0 -3556282.88386365 0\n"
]
}
],
"source": [
"!gdaltransform --help\n",
"# order to send in is Longitude Latitude Z\n",
"\n",
"# command-line gdaltransform -- lame command-line parsing but it works\n",
"metersXYZ = !echo 0.0 -60.0 | gdaltransform -s_srs \"$Mars_truncatedSphere\" -t_srs \"$srs\"\n",
"metersXYZ = str(metersXYZ).split(\"'\")[1]\n",
"X,Y,Z = str(metersXYZ).split()\n",
"\n",
"print(X,Y,Z)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['0 -60.0000000000001 0']\n"
]
}
],
"source": [
"# now go backwards to degrees\n",
"\n",
"#!gdaltransform --help\n",
"degreesLonLatZ = !echo $X $Y | gdaltransform -s_srs \"$srs\" -t_srs \"$Mars_truncatedSphere\" \n",
"print(degreesLonLatZ)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0 -3556282.8838636456 0.0\n"
]
}
],
"source": [
"# Now in Python, you can also use TransformPoint\n",
"\n",
"s_SRS = osr.SpatialReference() # makes an empty spatial ref object\n",
"s_SRS.ImportFromWkt(wkt_proj) # populates the spatial ref object with our WKT SRS\n",
"\n",
"t_SRS = osr.SpatialReference() # makes an empty spatial ref object\n",
"t_SRS.ImportFromProj4(Mars_truncatedSphere) # populates the spatial ref using Proj4 SRS\n",
"\n",
"# create the transform - this can be used repeatedly\n",
"transform = osr.CoordinateTransformation(t_SRS, s_SRS)\n",
"\n",
"# transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`\n",
" # order: Lon , Lat , Z\n",
"X,Y,Z = transform.TransformPoint(0.0, -60.0, 0.0)\n",
"\n",
"print(X,Y,Z)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['Report:',\n",
" ' Location: (1440P,1199L)',\n",
" ' Band 1:',\n",
" ' Value: 0.124235242605209']\n",
"\n",
"['Report:',\n",
" ' Location: (1440P,720L)',\n",
" ' Band 1:',\n",
" ' Value: 0.161960691213608']\n"
]
}
],
"source": [
"# Bonus section (for an X,Y, get pixel value)\n",
"# By using gdallocationinfo, see pixel location (row/col) and the value for that pixel\n",
"\n",
"import pprint\n",
"\n",
"pixel = !gdallocationinfo -geoloc $global_img $X $Y\n",
"pprint.pprint(pixel)\n",
"print()\n",
"\n",
"# try the center of the image at X,Y = 0,0\n",
"pixel = !gdallocationinfo -geoloc $global_img 0 0\n",
"pprint.pprint(pixel)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## crop to virtual file (VRT) during gdal_translate\n",
"\n",
"Now that we know the value in meters to crop with, we can crop out the south portion of the global image prior to re-projecting to south Polar Stereographic."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-10670400.0 5335200.0\n",
"10670400.0 -5335200.0\n",
"Input file size is 2880, 1440\n"
]
}
],
"source": [
"# now crop using meter bounding coordinates\n",
"out_crop = global_img + \"crop.vrt\"\n",
"\n",
"from osgeo import gdal\n",
"ds = gdal.Open(global_img)\n",
"gt = ds.GetGeoTransform() # captures origin and pixel size\n",
"\n",
"# upper left corner\n",
"ulx, uly = gdal.ApplyGeoTransform(gt,0,0)\n",
"print (ulx, uly)\n",
"\n",
"# lower right corner\n",
"lrx, lry = gdal.ApplyGeoTransform(gt,ds.RasterXSize,ds.RasterYSize)\n",
"print(lrx, lry)\n",
"ds = None\n",
"\n",
"# run gdal_translate to crop to a virtual file (*.vrt) on the command-line\n",
"# note we are using the $Y variable, or -60 degrees Latitude (in meters) from above\n",
"!gdal_translate -of vrt $global_img $out_crop -projwin $ulx $Y $lrx $lry"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_crop"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**After all that** -- a new option can do this in one command now during the gdaltranslate cropping stage. See the option \"-projwin_srs\" "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Input file size is 2880, 1440\n"
]
}
],
"source": [
"!gdal_translate -of vrt $global_img $out_crop -projwin -180 -60 180 -90 -projwin_srs \"$Mars_truncatedSphere\""
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_crop"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"original global image size is:\n",
"2880 columns x 1440 rows\n",
"\n",
"cropped image size is:\n",
"2880 columns x 240 rows\n",
"\n"
]
}
],
"source": [
"# Just for fun, using Python, look at the size difference between the original and cropped\n",
"# on the commandline - just use gdalinfo\n",
"\n",
"import json\n",
"import pprint\n",
"\n",
"# Open a original image dataset\n",
"dataset = gdal.Open(global_img, gdal.GA_ReadOnly)\n",
"\n",
"# How many rows and columns?\n",
"rows = dataset.RasterYSize\n",
"cols = dataset.RasterXSize\n",
"print('original global image size is:\\n{c} columns x {r} rows\\n'.format(c=cols, r=rows))\n",
"dataset = None\n",
"\n",
"# Open a crop image dataset\n",
"dataset = gdal.Open(out_crop, gdal.GA_ReadOnly)\n",
"\n",
"# How many rows and columns?\n",
"rows = dataset.RasterYSize\n",
"cols = dataset.RasterXSize\n",
"print('cropped image size is:\\n{c} columns x {r} rows\\n'.format(c=cols, r=rows))\n",
"dataset = None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## projecting to polar stereographic\n",
"\n",
"Why did we spend so much time cropping the image? This helps with extents for polar sterographic. We will try to use the global file and see the problems we have. Then we will try a cropped version.\n",
"\n",
"_gdalwarp_ is the application to do the map projection. https://gdal.org/programs/gdalwarp.html"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"+proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs \n",
"\n",
"Creating output file that is 2275P x 2279L.\n",
"Processing tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif.\n",
"Copying nodata values from source tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif to destination tes//out_south_pole.tif.\n",
"...10...20...30...40...50...60...70...80...90...100 - done.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n"
]
}
],
"source": [
"# get a south pole map projection from SpatialReference website\n",
"# note this site will return the recommended Mars radius (see the proj4 string below)\n",
"\n",
"import requests\n",
"response = requests.get('https://spatialreference.org/ref/iau2000/49920/proj4/')\n",
"IAU49920 = response.content.decode(\"utf-8\")\n",
"print(IAU49920)\n",
"print()\n",
"# now project the global file to polar\n",
"\n",
"#!gdalwarp\n",
"!gdalwarp -overwrite $global_img $out_SPole_img -t_srs \"$IAU49920\" -r bilinear"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_SPole_img"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**argh that image is horrible**. Now try to gdalwarp the virtual cropped version"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Creating output file that is 5801P x 5801L.\n",
"Processing tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tifcrop.vrt [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tifcrop.vrt.\n",
"Copying nodata values from source tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tifcrop.vrt to destination tes//out_south_pole.tif.\n",
"...10...20...30...40...50...60...70...80...90...100 - done.\n"
]
}
],
"source": [
"#!gdalwarp\n",
"!gdalwarp -overwrite $out_crop $out_SPole_img -t_srs \"$IAU49920\" -r bilinear"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_SPole_img"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Tip: if you ever see a gap in the polar image (or other projected image), you can use -wo SOURCE_EXTRA, for example -wo SOURCE_EXTRA=100. If the value used is too large you might see artifacts -- so watch out for that.\n",
"\n",
"example:\n",
"\n",
"> gdalwarp -overwrite in.tif out.tif -t_srs \"$IAU49920\" -r bilinear -wo SOURCE_EXTRA=100"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sinusoidal"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs \n",
"\n",
"Creating output file that is 2920P x 1357L."
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Processing tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif.\n",
"Copying nodata values from source tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif to destination tes//out_sinu_0.tif.\n",
"...10...20...30...40...50...60...70...80...90...100 - done.\n"
]
}
],
"source": [
"# Sinusoidal with a center at 0 Longitude\n",
"out_Sinu_img0 = \"tes//out_sinu_0.tif\"\n",
"\n",
"import requests\n",
"response = requests.get('https://spatialreference.org/ref/iau2000/49915/proj4/')\n",
"IAU49915 = response.content.decode(\"utf-8\")\n",
"print(IAU49915 + \"\\n\")\n",
"\n",
"#!gdalwarp\n",
"!gdalwarp -overwrite $global_img $out_Sinu_img0 -t_srs \"$IAU49915\" -r bilinear"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_Sinu_img0"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Creating output file that is 2920P x 1357L.\n",
"Processing tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif.\n",
"Copying nodata values from source tes//Mars_MGS_TES_Albedo_mosaic_global_7410m.tif to destination tes//out_sinu_180.tif.\n",
"...10...20...30...40...50...60...70...80...90...100 - done.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n",
"ERROR 1: Invalid dfSouthLatitudeDeg\n"
]
}
],
"source": [
"# Sinusoidal with a center at 180 Longitude\n",
"out_Sinu_img180 = \"tes//out_sinu_180.tif\"\n",
"\n",
"# update the definition above using a center meridian = 180 (+lon_0=180)\n",
"Sinu180 = \"+proj=sinu +lon_0=180 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs \"\n",
"\n",
"#!gdalwarp\n",
"!gdalwarp -overwrite $global_img $out_Sinu_img180 -t_srs \"$Sinu180\" -r bilinear"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"!tuiview $out_Sinu_img0 $out_Sinu_img180"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### bonus - test threading during gdalwarp, using HiRISE\n",
"\n",
"**Data for lesson**\n",
"\n",
"HiRISE file (copy in hirise// directory):\n",
" https://www.uahirise.org/PDS/DTM/ESP/ORB_045900_045999/ESP_045994_1985_ESP_046060_1985/ESP_046060_1985_RED_A_01_ORTHO.JP2\n",
" \n",
"recall for HiRISE we still may need to fix the GeoJpeg2000 header to correctly locate the image. For this timing exercise, we can skip the \"fix_jp2\" step:\n",
"https://planetarygis.blogspot.com/2016/07/more-hirise-conversion-tips-until.html"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"hirise = \"hirise//ESP_046060_1985_RED_A_01_ORTHO.JP2\"\n",
"to_hirise = \"hirise//out.tif\"\n",
"to_hirise2 = \"hirise//out2.tif\"\n",
"\n",
"# define Mars in degrees\n",
"Mars2000 = \"+proj=longlat +a=3396190 +b=3396190 +no_defs\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#Note here I am forcing 1 thread since gdalwarp will try to always be multi-threaded\n",
"\n",
"import time\n",
"\n",
"start = time.time()\n",
"!gdalwarp -overwrite $hirise $to_hirise -t_srs \"$Mars2000\" -r bilinear -wo NUM_THREADS=1\n",
"end = time.time()\n",
"print()\n",
"print(end - start)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# multithreaded using \"-wo NUM_THREADS=ALL_CPUS.\n",
"\n",
"start = time.time()\n",
"!gdalwarp -overwrite $hirise $to_hirise2 -t_srs \"$Mars2000\" -r bilinear -multi -wo NUM_THREADS=ALL_CPUS\n",
"end = time.time()\n",
"print()\n",
"print(end - start)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Speed summary. This may not be a huge speed-up for say a laptop, but on a more powerful machine (especially cluster machines with 24+ cpus), this setting can dramatically improve the speed."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment