Using Python gdalwarp/rasterio to read gif file and world file, re-project to another picture [closed]

Question

I want to re-project a radar image (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gif) which is in NAD83/EPSG4326 to WGS84/Pseudo-Mercator/EPSG3857. There is a world file associated with this gif (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gfw). I know I can run the following cmd command to re-project the picture:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 latest_radaronly.gif latest_radaronly_processed.tiff

I'm very new to GIS. I have two questions:

Question 1 Does the world file need to have exactly the same name as gif file?

Question 2 How can I do the re-reprojection using Python gdalwarp or rasterio package, and how to extract the coordinates bounds from the new re-projected picture?

I tested Marcelo Villa'answer and have two follow up questions:

Follow up question 1: How can I convert the coordinates in the bbox to longitudes (-180 to 180) and latitudes(-90 to 90)?

Follow up question 2: The tif generated is too large to be used in folium. I tried to resave the tif file to a gif file, however the map now looks like the following:

there is shade in the map. I was wondering what is the best way to reduce the picture size and avoid the shade in the map?

Answer 1

You can use the gdal.Warp() function to accomplish this. First save both files you linked to the same folder in your computer. Note that I had to create the .gfw file using Notepad ++ because when saving it from the URL it would save as a .txt file instead of a .gfw file. This is very important because otherwise gdal won't recognize the file and the raster won't be georeferenced. Then, you have to specify the output format for the reprojected file, the input spatial reference and the output spatial reference. If you want to specify other parameters, check the documentation. Finally, you just have to compute the bounding box of the reporjected image. This can be easily accomplished by taking the raster origin (left upper corner) and computing the right lower corner.

Here is a snippet to guide you:

import os

from osgeo import gdal


def compute_bbox(ds):
    """
    Computes a raster bounding box.

    Parameters
    ----------
    ds : gdal.Dataset object
        Raster to compute the bounding box for.

    Returns
    -------
    xmin : float
        Upper left x coordinate.
    ymin : float
        Lower right y coordinate.
    xmax : float
        Lower right x coordinate.
    ymax : float
        Upper left y coordinate.
    """
    width = ds.RasterXSize
    height = ds.RasterYSize
    xmin, pw, xskew, ymax, yskew, ph = ds.GetGeoTransform()
    xmax = xmin + (width * pw)
    ymin = ymax + (height * ph)

    return xmin, ymin, xmax, ymax


os.chdir(r'C:\path\to\folder')

kwargs = {'format': 'GTiff', 'srcSRS': 'EPSG:4326', 'dstSRS': 'EPSG:3857'}
ds = gdal.Warp('latest_radaronly_3857.tif', 'latest_radaronly.gif', **kwargs)

bbox = compute_bbox(ds)

If you inspect bbox you will get:

>>>bbox
(-14207635.496435506, 2470074.029222458, -7406084.60511778, 6518566.074162895)

Don't forget to close the dataset at the end of the script with the following line:

ds = None

---

Answering to your follow-up questions.

1) Reprojecting the image to EPSG 3857 results in a file that uses meters instead of degrees as its unit. If you want to get the bounding box coordinates in degrees rather than meters, you can execute the get_bbox() function on the original dataset like this:

from osgeo import gdal

ds = gdal.Open('latest_radaronly.gif', 0)
bbox = get_bbox(ds)

2) You can reduce the size of the resulting TIFF file by specifying a type of compression. Take a look at several questions on this site regarding that matter.

On another note, avoid asking multiple questions in one single question as well as editing your original question to ask further questions.