I have created tif image using matplotlib python trigulation function, Now tif need to be change to goetiff with epsg4326 projection with removing the background only the patch need to be there.
What difficulty i am facing :- 1- projection change is not happening correct. 2- pixel value should be z parameter value. 3- only the clip of lat lon i want not the whole lat lon area.
CODES:- 1- For creating matplotlib trigulation tiff code below
#MAKING TIFF FILE USING TRIANGULARION MATPLOTLIB.
import matplotlib.tri as mtri
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def plot():
x=np.random.uniform(70,94,74)
y=np.random.uniform(8,22,74)
# A selected triangulation of the points.
triangles = np.asarray([
[67, 66, 1], [65, 2, 66], [1, 66, 2], [64, 2, 65], [63, 3, 64],
[60, 59, 57], [2, 64, 3], [3, 63, 4], [0, 67, 1], [62, 4, 63],
[57, 59, 56], [59, 58, 56], [61, 60, 69], [57, 69, 60], [4, 62, 68],
[6, 5, 9], [61, 68, 62], [69, 68, 61], [9, 5, 70], [6, 8, 7],
[4, 70, 5], [8, 6, 9], [56, 69, 57], [69, 56, 52], [70, 10, 9],
[54, 53, 55], [56, 55, 53], [68, 70, 4], [52, 56, 53], [11, 10, 12],
[69, 71, 68], [68, 13, 70], [10, 70, 13], [51, 50, 52], [13, 68, 71],
[52, 71, 69], [12, 10, 13], [71, 52, 50], [71, 14, 13], [50, 49, 71],
[49, 48, 71], [14, 16, 15], [14, 71, 48], [17, 19, 18], [17, 20, 19],
[48, 16, 14], [48, 47, 16], [47, 46, 16], [16, 46, 45], [23, 22, 24],
[21, 24, 22], [17, 16, 45], [20, 17, 45], [21, 25, 24], [27, 26, 28],
[20, 72, 21], [25, 21, 72], [45, 72, 20], [25, 28, 26], [44, 73, 45],
[72, 45, 73], [28, 25, 29], [29, 25, 31], [43, 73, 44], [73, 43, 40],
[72, 73, 39], [72, 31, 25], [42, 40, 43], [31, 30, 29], [39, 73, 40],
[42, 41, 40], [72, 33, 31], [32, 31, 33], [39, 38, 72], [33, 72, 38],
[33, 38, 34], [37, 35, 38], [34, 38, 35], [35, 37, 36]])
z = np.random.uniform(0.5, 8, 74)
triang = mtri.Triangulation(x, y, triangles=triangles)
print("type(triang):", triang)
print("type(z):", type(z))
ax=plt.tripcolor(triang,z,vmin=0,vmax=2)
ax
plt.savefig("ax.tif")
plt.show(ax)
plot()
2- making geotiff with projection from tif
from osgeo import gdal, osr
import numpy, gdal
ds = gdal.Open('ax.tif')
band = ds.GetRasterBand(1)
myarray = numpy.array(band.ReadAsArray())
convert=lambda x:x*(0.0124796)
myarray=convert(myarray)
#print("myarray:",myarray)
selection = numpy.logical_or(myarray >= 3.5, myarray <= 7)
new_array = [ [0 for i in range(band.XSize)] for j in range(band.YSize)]
for i, item in enumerate(myarray):
for j, element in enumerate(item):
if selection[i][j] == True:
new_array[i][j] = myarray[i][j]
else:
new_array[i][j] = -999
if new_array[i][j]==3.182298:
new_array[i][j] = -999
geotransform = ds.GetGeoTransform()
wkt = ds.GetProjection()
# Create gtif file
driver = gdal.GetDriverByName("GTiff")
output_file = "ax_geotiff.tif"
dst_ds = driver.Create(output_file,
band.XSize,
band.YSize,
1,
gdal.GDT_Float64)
new_array = numpy.array(new_array)
print(type(new_array[i]))
#new_array[i]=new_array.max())
#transformed_arr = np.where(arr<5, arr*10, 0)
#new_array=numpy.where(new_array[i]==new_array.max(),-999,new_array)
print("new_array:",new_array)
#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(-999)
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
#srs = osr.SpatialReference()
#srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset
ds = None
dst_ds = None
print(numpy.unique(new_array))
3- To plot the geotiff
#CHECKING THE TIFF FILE / ANALYSING THE TIFF FILE DATA AND PLOT.
#tttraster_geotiff
import gdal
%matplotlib inline
import rasterio
from matplotlib import pyplot
import numpy
ds = gdal.Open('ax_geotiff.tif')
band = ds.GetRasterBand(1)
myarray = numpy.array(band.ReadAsArray())
print(myarray)
pyplot.imshow(myarray, cmap='pink')
pyplot.show()
##Above is the sample data and codes i am using but i am not getting correct epsg 4326 projection and my background color in geottiff blue color is coming if i select cmap='jet' i.e blue color large part of geotiff plot indicate the mav value of z parameter , those blue color if i change to np.nan then i a not able to see the other,
so what all modification i have to do- 1- correct projection 2- remove only the blue color as nan color.enter image description here enter image description here 3- want to remove that drak red color(max value of z variable) , only the map plot i want.
