1

I want to reproject a raster using the projection of other raster using gdalwarp, however, I am confronting some problems because the raster does not line up when loading them into a GIS Software.

This is my code:

from osgeo import gdal,osr
out_ras_500m_buff = r'MK_020m_nat_buffer500m_v2.tif'
boundaries_in = r'MK_020m_nat_buffer100m_v2.tif'
ds=gdal.Open(boundaries_in)
prj=ds.GetProjection()
print(prj)
gt = ds.GetGeoTransform()
cols = ds.RasterXSize 
rows = ds.RasterYSize

srs=osr.SpatialReference(wkt=prj)
if srs.IsProjected:
    print(srs.GetAttrValue('projcs'))
print(srs.GetAttrValue('geogcs'))

ext = GetExtent(gt, cols, rows)
print(ext)
xmin = ext[0][0]
ymin = ext[2][1]
xmax = ext[2][0]
ymax = ext[0][1]
print(f"XMin:{xmin}")
print(f"YMin:{ymin}")
print(f"XMax:{xmax}")
print(f"YMax:{ymax}")

EPSG_CODE = "6204"
warped_raster = out_ras_500m_buff.replace(".tif","_" + EPSG_CODE+ ".tif")
warped_raster = warped_raster.replace("02", "03")

prj_file = OutDir + "/" + os.path.basename(out_shp).replace(".shp", ".wkt")

with open(prj_file, "w") as f:
    f.write(prj)
#subprocess.call("gdalsrsinfo -o wkt {raster_boundaries_nat} > target.wkt")

#gdal.Warp(warped_raster, out_ras_500m_buff, dstSRS= f"EPSG:{EPSG_CODE}")

#subprocess.call(f'gdalwarp -t_srs {prj} {out_ras_500m_buff} {warped_raster} -te {xmin} {ymin} {xmax} {ymax}')
command =  f'gdalwarp -t_srs {prj_file} {out_ras_500m_buff} {warped_raster} -te {xmin} {ymin} {xmax} {ymax} -tr 10 10 -r near -co "COMPRESS=LZW"'
subprocess.call(command)

I am using a GetExtent function which get the extent of the raster I want to use to reproject the other one:

def GetExtent(gt,cols,rows):
    ''' Return list of corner coordinates from a geotransform

        @type gt:   C{tuple/list}
        @param gt: geotransform
        @type cols:   C{int}
        @param cols: number of columns in the dataset
        @type rows:   C{int}
        @param rows: number of rows in the dataset
        @rtype:    C{[float,...,float]}
        @return:   coordinates of each corner
    '''
    ext=[]
    xarr=[0,cols]
    yarr=[0,rows]

    for px in xarr:
        for py in yarr:
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
        yarr.reverse()
    return ext

Although the result contains the CRS that I want, the projected raster does not have the following information that is found in a .prj file for the target spatial reference:

PROJCS["Macedonian_State_Coordinate_System",
    GEOGCS["GCS_MGI_1901",
        DATUM["D_MGI_1901",
            SPHEROID["Bessel_1841",6377397.155,299.1528128],
            TOWGS84[-521.7476,-229.4892,-509.9207,-4.02878,-4.48836,15.52067,-9.7803]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",21.0],
    PARAMETER["Scale_Factor",0.9999],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0],
    AUTHORITY["EPSG",6204]]

My question is:

TOWGS84[-521.7476,-229.4892,-509.9207,-4.02878,-4.48836,15.52067,-9.7803]]

Given that information, how can I use the .prj file to reproject a raster with the correct parameters (shift to WGS84)?

It seems that only giving the EPSG CODE in the subprocess.call does not make the work.

I found that using pyproj, and setting a string like the following, one can accomplish that. However, I cannot understand how to plug the shift to wg84 parameters when reprojection a raster.

isn2004=pyproj.Proj("+proj=lcc +lat_1=64.25 +lat_2=65.75 +lat_0=65 +lon_0=-19 +x_0=1700000 +y_0=300000 +no_defs +a=6378137 +rf=298.257222101 +to_meter=1")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.