0

I'm new in QGIS and I’m having some troubles. I have produced a surface plot in python and I would like to import it in QGIS but I don't know how do that, can anyone help me?

P.S. I also have the long lat coordinates of the plot for the geo referencing

data = pandas.read_csv('TEST.csv', sep = ';', header = None)
data = numpy.array(data,dtype = numpy.float)

AZI = data[:,0]
LAT = data[:,1]
LON = data[:,2]
DIST = data[:,3]
EdB = data[:,4]
HBS_STEP = data[:,5]

TxCoor = data[0, 1:3];

minHBS = int(numpy.min(HBS_STEP))
maxHBS = int(numpy.max(HBS_STEP))

DIST_NEW = numpy.zeros([(maxHBS + 1), 360])
EdB_NEW = numpy.zeros([(maxHBS + 1), 360])

ELE_NEW = numpy.zeros([(maxHBS + 1), 360])

AZI_V = numpy.linspace(0, 359, 360)
AZI_NEW = numpy.matlib.repmat(AZI_V, 36, 1)

indexVector = numpy.empty(0)
distMax = 0

count = 0

rangeMinMax = range(minHBS, (maxHBS + 1))

for i in rangeMinMax:

    indexVector = numpy.empty(0)
    distMax = 0

    index = numpy.where(HBS_STEP == i)
    AZIs = AZI[index]

    for j in range(0, len(AZIs)):

        if (j == 0 and index[0][j] != index[0][j+1]):

            DIST_NEW[i][int(AZIs[j])] = DIST[index[0][j]];
            EdB_NEW[i][int(AZIs[j])] = EdB[index[0][j]];

        elif (j > 0 and index[0][j-1] != index[0][j] and len(indexVector) == 0):

            DIST_NEW[i][int(AZIs[j])] = DIST[index[0][j]]
            EdB_NEW[i][int(AZIs[j])] = EdB[index[0][j]]

        elif (j > 0 and index[0][j-1] == index[0][j]):

            indexVector = numpy.concatenate(indexVector, index)

        else:

            count = count + 1

            for z in range(0, len(indexVector)):

                if(DIST[int(AZIs[j])] > distMax):

                    distMax = DIST[AZIs[j]]

                    DIST_NEW[i][int(AZIs[j])] = DIST[index[0][j]]
                    EdB_NEW[i][int(AZIs[j])] = EdB[index[0][j]] 

            indexVector = numpy.empty(0)
            distMax = 0

for w in range(0, len(DIST_NEW[:,1])):
    for v in range(0, len(DIST_NEW[1,:])):
        if (DIST_NEW[w, v] == 0):
            ELE_NEW[w, v] = numpy.nan
        else:
            ELE_NEW[w, v] = numpy.round(numpy.arctan((w-1)/(DIST_NEW[w,v]*100))*(180/numpy.pi)*10)/10
#            numpy.round((((numpy.arctan((w-1)/(numpy.multiply(numpy.multiply((DIST_NEW(w,v),1000)),180))))/numpy.pi)*10))/10

numpy.savetxt('ELE_NEW.txt', ELE_NEW, fmt="%1.2e", delimiter=',')
numpy.savetxt('AZI_NEW.txt', AZI_NEW, fmt="%1.2e", delimiter=',')
numpy.savetxt('DIST_NEW.txt', DIST_NEW, fmt="%1.2e", delimiter=',')

ELE_NEW = (numpy.pi/2)+(ELE_NEW/180)*(numpy.pi)
AZI_NEW = (AZI_NEW/180)*(numpy.pi)

xconv = DIST_NEW*numpy.sin(ELE_NEW)*numpy.cos(AZI_NEW)
yconv = DIST_NEW*numpy.sin(ELE_NEW)*numpy.sin(AZI_NEW)
zconv = DIST_NEW*numpy.cos(ELE_NEW)

############################################################################

fig = plt.figure()

# ax = fig.add_subplot(1,1,1, projection='3d') 
ax = fig.gca(projection='3d')
ax.grid(True)
ax.axis('on')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# colormap
minn, maxx = EdB_NEW.min(), EdB_NEW.max()
norm = matplotlib.colors.Normalize(minn, maxx)
m = plt.cm.ScalarMappable(norm=norm, cmap='jet')
m.set_array([])
fcolors = m.to_rgba(EdB_NEW)     
surf = ax.plot_surface(xconv, yconv, zconv, rstride=1, cstride=1, facecolors=fcolors, vmin=minn, vmax=maxx)

plt.show()

plt.close(fig)


enter image description here

  • You need to add more details. Your data is in what format? Numpy Array? Add a reproducable code snippet. – BERA Feb 14 at 9:22
  • Hi BERA, thank you for the answer. I have added the code and the result but i can't upload the csv file. Anyway I also have the long lat coordinates of every point. – FM79 Feb 14 at 9:43
  • @BERA maybe i can add xconv, yconv and zconv values if you think that may be usefull – FM79 Feb 14 at 11:24
  • Do you want to turn it into raster or vector data in QGIS? What is your next step – BERA Feb 14 at 11:43
  • @BERA Maybe a raster data is better – FM79 Feb 14 at 11:49
1

Using Writing numpy array to raster file as suggested by BERA is the easiest way

With one of my examples

Matplotlib 2D:

enter image description here

Matplotlib 3D (with def axisEqual3D(ax) in set matplotlib 3d plot aspect ratio):

ax.plot_surface(xconv, yconv, zconv, rstride=1, cstride=1, cmap='gist_earth',antialiased=True)
ax.view_init(60,-160)
axisEqual3D(ax)

enter image description here

Convert it into an ASCII Raster

xmin,ymin,xmax,ymax = [xconv,.min(),yconv.min(),xconv.max(),yconv.max()]
nrows,ncols = np.shape(zconv)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)  
# creation of the raster file 
header = "ncols     %s\n" % ncols
header += "nrows    %s\n" % nrows
header += "xllcorner %s\n" % geotransform[0]
header += "yllcorner %s\n" % geotransform[3]
header += "dx  %s\n" % geotransform[1]
header += "dy  %s\n" % geotransform[5]
header += "NODATA_value -9999\n"
with  open("grid.asc", "w") as f:
  f.write(header)
  np.savetxt(f,zconv, fmt="%1.2f")

grid.asc content

ncols 100
nrows 100
xllcorner 251215.52097430476
yllcorner 46722.44250503713
dx 11.230954009630368
dy -8.647123208649791
NODATA_value -9999
395.37 395.23 395.10 ...

Result in QGIS:

enter image description here

But you can also use osgeo.gdal (Writing Numpy array to raster file (tif) returns a trivial black square) or rasterio (Lesson 1. Export Numpy Arrays to Geotiff Format Using Rasterio and Python)

|improve this answer|||||

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.