# Converting coordinates to pixels with out losing points

I'm trying to convert coordinates from lat/lon into pixel but I lose points in this process.

The code that I'm using is the following:

``````from matplotlib import pyplot as plt
import numpy as np
import pandas as pd

cords=df.as_matrix(columns=['x','y'])
gt=[7.6445503225, 5.4065168747250134e-06,  0.0,  45.07436634583334,  0.0, -5.406516856707135e-06]

index=np.zeros(cords.shape)
index[:,1]=((cords[:,1] - gt) / gt).round()
index[:,0]=((cords[:,0] - gt) / gt).round()
index=index.astype(int)
index[:,0]=index[:,0]-min(index[:,0])+1
index[:,1]=index[:,1]-min(index[:,1])+1
row=max(index[:,1])
col=max(index[:,0])
image=np.zeros([row+1,col+1])
for i in range(0,len(index)):
image[index[i,1],index[i,0]]=255
``````

If I plot the cords or the index points I get this: If I plot the image I get this: As you can see there are some points that are missing in translating lat/lon in to pixel numbers. Yellow is 255 value and purple is 0 value. How can this be solved?

Here you find the coordinates that I'm using cords.csv

Here you find the coordinates with the values that need to be set to each pixel. cords_valus.csv

EDIT

In the folowing link you find the raster image from where I have extract the cordinates. It fits to Italy perfectly. raster

• link to coordinate file is broken – Stephen Rauch Jul 29 '17 at 21:54
• I think that there is not any problem in your data. They have values equal zero in that area. Issues are others. Please, see my answer. – xunilk Jul 29 '17 at 23:19

I think that there is not any problem in your data. They have values equal zero in that area. Issues are others. When I modified your code for obtaining your points and your raster (by using gdal and PyQGIS python modules), they don't match in Italy area. Your geotransform parameters are wrong (I fixed them with 'Polygon from layer extent' for point layer). Complete code is as follow:

``````from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from osgeo import gdal, osr

cords = df.as_matrix(columns=['x','y'])

points = [ QgsPoint(pt, pt) for pt in cords ]

#gt = [7.6445503225, 5.4065168747250134e-06,  0.0,  45.07436634583334,  0.0, -5.406516856707135e-06]
gt = [7.65917, 5.4065168747250134e-06,  0.0,  45.0659,  0.0, -5.406516856707135e-06]
#xMin,yMin 7.65917,45.0656 : xMax,yMax 7.65972,45.0659

index=np.zeros(cords.shape)

index[:,1]=((cords[:,1] - gt) / gt).round()
index[:,0]=((cords[:,0] - gt) / gt).round()

index=index.astype(int)

index[:,0]=index[:,0]-min(index[:,0])+1
index[:,1]=index[:,1]-min(index[:,1])+1

row=max(index[:,1])
col=max(index[:,0])

image=np.zeros([row+1,col+1])

for i in range(0,len(index)):
image[index[i,1],index[i,0]] = df['value'][i]

# Create gtif file
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/image.tif"

dst_ds = driver.Create(output_file,
col+1,
row+1,
1,
gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( image )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(gt)

dst_ds.GetRasterBand(1).SetNoDataValue(0)

# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
dst_ds.SetProjection( srs.ExportToWkt() )

dst_ds = None

epsg = 4326

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:double""&index=yes"

mem_layer = QgsVectorLayer(uri,
'point',
'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
feat.setAttributes([ i, float(df['value'][i]) ])
feat.setGeometry(QgsGeometry.fromPoint(points[i]))

``````

After running above code at Python Console, it can be observed that points and image are shifted. It's due to algorithm to produce image. You need to fix it. However, you can also rasterize your points by values and get a very quick result; as it observed at next image: Apparently, there is not lost points in that area; as it showed by above image.

• I have obtained the gt vectort from gdal.geotrasform form the raster data. the points in cords.csv are obtained from a raster image that don't have that line in the middle. The raster image (yellow) is obtained from the points after chanign the coordinates from wgs84 to pixel coordinates. I belive that the problem is ralated to the round(). – Lorenzo Bottaccioli Jul 29 '17 at 23:25
• Check the edit for the raster image. – Lorenzo Bottaccioli Jul 30 '17 at 0:04