Friday, February 14, 2025

AutoLISP Program that iterates through selected lines, calculates their angles and lengths, and then annotates them

 In other words, the autoLISP program will select a line, get its bearing and distance and write the values on beside the line. The distance is in meters and the bearing format is in Whole Circle Bearing and DD°MM'SS


(defun c:BRGDIST (/ ss i entData ptStart ptEnd dx dy angleRad angleDeg wcb bearingStr dist midpoint)
  ; Function to convert decimal degrees to DMS (degrees, minutes, seconds)
  (defun DEGtoDMS (deg / d m s)
    (setq d (fix deg)) ; Degrees
    (setq m (fix (* (- deg d) 60))) ; Minutes
    (setq s (* (- (* (- deg d) 60) m) 60)) ; Seconds
    (strcat (itoa d) "°" (itoa m) "'" (rtos s 2 2) "\"") ; Format as DD°MM'SS"
  )

  ; Select lines
  (setq ss (ssget '((0 . "LINE"))))
  (if ss
    (progn
      (setq i 0)
      (repeat (sslength ss)
        ; Extract line data
        (setq entData (entget (ssname ss i)))
        (setq ptStart (cdr (assoc 10 entData))) ; Start point
        (setq ptEnd (cdr (assoc 11 entData)))   ; End point
        
        ; Calculate bearing
        (setq angleRad (angle ptStart ptEnd)) ; Angle in radians (from X-axis, counterclockwise)
        (setq angleDeg (/ (* angleRad 180.0) pi)) ; Convert to degrees
        
        ; Convert to Whole Circle Bearing (clockwise from north, 0°-360°)
        (setq wcb (rem (+ (- 90 angleDeg) 360) 360)) ; Adjust for WCB
        (setq bearingStr (DEGtoDMS wcb)) ; Convert WCB to DMS format
        
        ; Calculate distance
        (setq dist (rtos (distance ptStart ptEnd) 2 2)) ; Distance with 2 decimal places
        
        ; Calculate midpoint
        (setq midpoint (list
          (/ (+ (car ptStart) (car ptEnd)) 2.0)
          (/ (+ (cadr ptStart) (cadr ptEnd)) 2.0)
        ))
        
        ; Create text at midpoint
        (entmake
          (list
            '(0 . "TEXT")
            (cons 10 midpoint)          ; Insertion point
            (cons 40 5)               ; Text height
            (cons 1 (strcat bearingStr "      " dist"m")) ; Bearing (DMS) + distance
            (cons 50 angleRad)          ; Rotation angle (aligned with line)
            (cons 7 "Standard")         ; Text style
          )
        )
        (setq i (1+ i))
      )
    )
  )
  (princ)
)


This will create a new layer for the annotation

(defun c:BRGDIST (/ ss i entData ptStart ptEnd dx dy angleRad angleDeg wcb bearingStr dist midpoint)
  ; Function to convert decimal degrees to DMS (degrees, minutes, seconds)
  (defun DEGtoDMS (deg / d m s)
    (setq d (fix deg)) ; Degrees
    (setq m (fix (* (- deg d) 60))) ; Minutes
    (strcat (itoa d) "°" (itoa m) "'") ; Format as DD° MM'
    ;; (setq s (* (- (* (- deg d) 60) m) 60)) ; Seconds
    ;; (strcat (itoa d) "°" (itoa m) "'" (rtos s 2 2) "\"") ; Format as DD° MM' SS"
  )

  ; Create or use the annotation layer
  (if (not (tblsearch "LAYER" "BRG_DIST"))
    (command "._LAYER" "_M" "BRG_DIST" "_C" "1" "" "") ; Create layer and set color to red
  )
  (setvar "CLAYER" "BRG_DIST") ; Set current layer to BRG_DIST

  ; Select lines
  (setq ss (ssget '((0 . "LINE"))))
  (if ss
    (progn
      (setq i 0)
      (repeat (sslength ss)
        ; Extract line data
        (setq entData (entget (ssname ss i)))
        (setq ptStart (cdr (assoc 10 entData))) ; Start point
        (setq ptEnd (cdr (assoc 11 entData)))   ; End point
        
        ; Calculate bearing
        (setq angleRad (angle ptStart ptEnd)) ; Angle in radians (from X-axis, counterclockwise)
        (setq angleDeg (/ (* angleRad 180.0) pi)) ; Convert to degrees
        
        ; Convert to Whole Circle Bearing (clockwise from north, 0°-360°)
        (setq wcb (rem (+ (- 90 angleDeg) 360) 360)) ; Adjust for WCB
        (setq bearingStr (DEGtoDMS wcb)) ; Convert WCB to DMS format
        
        ; Calculate distance
        (setq dist (rtos (distance ptStart ptEnd) 2 2)) ; Distance with 2 decimal places
        
        ; Calculate midpoint
        (setq midpoint (list
          (/ (+ (car ptStart) (car ptEnd)) 2.0)
          (/ (+ (cadr ptStart) (cadr ptEnd)) 2.0)
        ))
        
        ; Create text at midpoint
        (entmake
          (list
            '(0 . "TEXT")
            (cons 10 midpoint)          ; Insertion point
            (cons 40 5)               ; Text height
            (cons 1 (strcat bearingStr "   " dist "m")) ; Bearing (DMS) + distance
            (cons 50 angleRad)          ; Rotation angle (aligned with line)
            (cons 7 "Standard")         ; Text style
            (cons 8 "BRG_DIST")         ; Layer
          )
        )
        (setq i (1+ i))
      )
    )
  )
  (princ)
)

Happy Coding!

Thursday, February 6, 2025

Generating AutoCAD Script (.SCR) file from GPS/GNSS Data

 Recently, I was working on a project where I had to label some points in several AutoCAD files from data points obatined from GPS/GNSS reciever. So, I wrote this Python code to 'Generate AutoCAD Script (.SCR) file from GNSS Data'.

easting = '356662'
northing = '946418'

# construct text string eg: _TEXT easting,northing 10 90 PA10
script_str_txt = f'_TEXT {easting},{northing} 10 90 {name}\n'

# construct point string eg: _POINT easting,northing
script_str_pt = f'_POINT {easting},{northing}\n'

with open(f"{file_num}_{sch_name}.scr", "a") as f:
    f.write(script_str_txt)

The result of the script:


A more practical application code is this one below. It reads a large CSV file containing school names, then groups the data by the school names and generate .scr file for each group.
sch_dict = {'GSS TOTO WEST TOTO LGC':1, 'GJSS ASOPADA KARU LGC':2, 'GJSS KURMIN TAGWAYE AKWANG':3, 'GJSS LAFIA SOUTH LAFIA LGC':4, 'GJSS SAMBGOBARAU KOKONA LGC':5, 'GOV DAY SEC OBI LGC':6, 'GS LOKO NASARAWA LGC':7, 'GSS ALUSHI KEANA LGC':8, 'GSS AZUB CENTRAL LAFIA LGC':9, 'GSS EFUGBORINGO DOMA LGC':10, 'GSS GUDI AKWANGA LGC':11, 'GSS KEKURA AWE LGC':12, 'GSS MAMA WAMBA LGC':13, 'GSS SABONGARI KEFFI KEFFI LGC':14, 'ISLAMIYYA LGEA KEANA LGC':15, 'LGEA DHIZILA OBI LGC':16, 'LGEA EFURIGBORINGO DOMA LGC':17, 'LGEA EZZEN MADA STATION EGGON':18, 'LGEA GUDIGE TSOHO NASARAWA LGC':19, 'LGEA IBI AWE LGC':20, 'LGEA KIGBUNA LAFIA LGC':21, 'LGEA MADA HILL AKWANGA LGC':22, 'LGEA PILOT CEN UMAISHA TOTO LG':23, 'LGEA PRI ADOKASA KARU LGC':24, 'LGEA PRI KONZA WAMB LGC':25, 'LGEA PRI NASARAWA EAST NAS LGC':26, 'LGEA REFUGE CAMP LAFIA LGC':27, 'LGEA SAKWATO KOKONA LGC':28, 'LGEA TRANSFER UMAISHA TOT LGC':29, 'NADP PROJECT LAFIA LGC':30, 'RCM GSS NAS EGGON LGC':31}


# Read the GNSS created spread sheet....
df = pd.read_csv(r"sch-24.csv")

group_df = df.groupby('School Name')
group_keys = list(group_df.groups.keys())

for s in group_keys:
    temp_df = group_df.get_group(s)

    for idx, row in temp_df.iterrows():
        easting = row['X']
        northing = row['Y']
        name = row['Name']
        sch_name = row['School Name']
        file_num = sch_dict[sch_name]
    
        # construct text string eg: _TEXT easting,northing 10 90 PA10
        script_str_txt = f'_TEXT {easting},{northing} 10 90 {name}\n'
        # construct point string eg: _POINT easting,northing
        script_str_pt = f'_POINT {easting},{northing}\n'

        with open(f"{file_num}_{sch_name}.scr", "a") as f:
            f.write(script_str_txt)

        with open(f"{file_num}_{sch_name}.scr", "a") as f:
            f.write(script_str_pt)
        
    # break

print('Done...')



Happy coding!

Thursday, January 30, 2025

Images Georeferencing World Files

 World Files are georeferencing files for images that describe the location, scale, and orientation of an image such geoTiff with .tfw as its world file, JPEG with .jgw as its world file and PNG with .pgw as its world file.

Image world files (.tfw, .jgw, .pgw) are text files that contains geographic information for a raster image. They are used to georeference the image, which means to position the image correctly on a map or in a mapping system.

The .tfw world file is a text file used to georeference the GeoTIFF raster images, like the orthomosaic and the DSM. The .jgw world file is a text file used to georeference the JPEG raster images. The .pgw world file is a text file used to georeference the PNG raster images.

Any of the world files have similar structure which contain 6-line in the text file:

  • Line 1: pixel size in the x-direction in map units (GSD).
  • Line 2: rotation about y-axis.
  • Line 3: rotation about x-axis.
  • Line 4: pixel size in the y-direction in map in map units (GSD).
  • Line 5: x-coordinate of the upper left corner of the image.
  • Line 6: y-coordinate of the upper left corner of the image.

*GSD: Ground Sampling Distance.

As pixels are considered as square lines 1 and 4 are the same.


Georeferencing a JPEG image in AutoCAD using its .jgw World Files

This one used case of the worl files in AutoCAD. The Visual LISP program below adopted from Insert georeferenced image with world file in AutoCAD will select an image that has a jgw world file and georeferencing it appropraitely.

(vl-load-com)
 (defun ss->lst (ss / e n out)
   (setq n -1)
   (while (setq e (ssname ss (setq n (1+ n)))) (setq out (cons (vlax-ename->vla-object e) out)))
 )
 (defun _writefile (filename lst / file result)
   (cond ((and (eq 'str (type filename)) (setq file (open filename "w")))
   (foreach x lst
     (write-line
       (cond ((= (type x) 'str) x)
	     ((= (type x) 'int) (itoa x))
	     ((= (type x) 'real) (rtos x 2 6))
	     ((vl-prin1-to-string x))
       )
       file
     )
   )
   (close file)
   filename
  )
   )
 )
 (defun _readfile (filename / file result)
   (cond
     ((and (eq 'str (type filename)) (setq file (open filename "r")))
      (while (setq line (read-line file)) (setq result (cons (vl-string-trim " " line) result)))
      (close file)
      (reverse result)
     )
   )
 )
 (setq opt "ReadIt")
;  (initget 0 "ReadIt WriteIt")
;  (setq	opt (cond ((getkword (strcat "\nImage World File [ReadIt/WriteIt] <" opt ">: ")))
;		  (opt)
;	    )
; )
 (princ "\nSelect image(s): ")
 (setq pre (getvar 'dwgprefix))
 (if (and (setq ss (ssget '((0 . "image")))) (setq ss (ss->lst ss)))
   (foreach image ss
     (setq name    (vlax-get image 'name)
    hgt	    (vlax-get image 'height)
    wdth    (vlax-get image 'width)
    imhgt   (vlax-get image 'imageheight)
    imwdth  (vlax-get image 'imagewidth)
    rot	    (vlax-get image 'rotation)
    bpt	    (vlax-get image 'origin)
    imgpath (vl-filename-directory (vlax-get image 'imagefile))
    jgw	    (strcat imgpath "\\" name ".jgw")
     )
     (if (= opt "ReadIt")
(progn
  (if (and (or (setq jgw (findfile (strcat pre name ".jgw")))
	       (setq jgw (findfile (strcat imgpath "\\" name ".jgw")))
	       (setq jgw (getfiled (strcat "***Select <<" name ".jgw>>***") pre "jgw" 16))
	   )
	   (setq pre (strcat (vl-filename-directory jgw) "\\"))
	   (setq data (mapcar 'atof (_readfile jgw)))
	   (> (length data) 5)
	   (setq l1 (car data))
	   (setq mvpt (list (nth 4 data) (nth 5 data) 0.0))
      )
    (progn (vla-put-imageheight image (* hgt l1))
	   (vla-put-imagewidth image (* wdth l1))
	   (vla-put-rotation image (cadr data))
	   (setq rot (vlax-get image 'rotation))
	   (setq bpt (polar bpt (+ (/ pi 2.) rot) (* hgt l1)))
	   (vlax-invoke image 'move bpt mvpt)
	   (princ (strcat "\njgw File Read - " jgw))
    )
    (princ "\njgw file NOT found or not correctly formatted!")
  )
)
(progn (setq bpt (polar bpt (+ (/ pi 2.) rot) imhgt))
       (if (setq jgw (_writefile
		       (strcat imgpath "\\" name ".jgw")
		       (list (/ imhgt hgt)
			     rot
			     (strcat "-" (rtos (abs rot) 2 6))
			     (strcat "-" (rtos (abs (/ imwdth wdth)) 2 6))
			     (rtos (car bpt) 2 6)
			     (rtos (cadr bpt) 2 6)
		       )
		     )
	   )
	(print jgw)
	(princ "\nError writing file...")
       )
)
     )
   )
 )
 (princ)

Save the script in a lisp file or just copy and paste it into the vlisp window as shown below. Then run the program and select the image you want to georeference in the AutoCAD current window. this will move, scale, rotate and orient the image to its real-world coordinates position. Zoom to Extend to see the result.


Thank you for reading.

Saturday, January 25, 2025

Stop QGIS from displaying attribute dialog box during digitizing

 In QGIS the attribute entry form pop-up immediately you digitized a point or line or polygon feature. This is so that you can enter the corresponding attribute to the feature immediately, that makes sense. Howerver, it is not all the time that you may want this to happend. There are time you just don't care about the attributes and all you require is the geometrical features didgitized asap. In such a situation the attribute dialog box pop-up during digitizing may be anoying.

Here is how to prevent the attribute dialog box from displaying during digitizing. You just need to enable the feature: "Suppress attribute from pop-up after feature creation".

Go to: Settings menu >> Options >> Map Tools >> Digitizing



That is it!

Saturday, January 4, 2025

GEE - Measuring Land Surface Temprature from Landsat 8

 The Javascript code below will calculate Land Surface Temperature (LST) from landsat 8 imagery.



Adopted from this medium post: Analyzing Land Surface Temperature (LST) with Landsat 8 Data in Google Earth Engine - https://medium.com/@ridhomuh002/analyzing-land-surface-temperature-lst-with-landsat-8-data-in-google-earth-engine-f4dd7ca28e70
// Load AOI or ROI
var aoi = ee.FeatureCollection("users/forum2k9/Karu"); // Karu polygon - change if you need different 
Map.addLayer(aoi, {}, 'AOI - Karu');
Map.centerObject(aoi, 10);

// Define data retrieval Start and End Date
var startDate = '2023-03-01';
var endDate = '2023-10-31';


// Applies scaling factors. Recalculating from K to °C
function applyScaleFactors(image) {
 // Scale and offset values for optical bands
 var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
 // Scale and offset values for thermal bands
 var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0).subtract(273.15); // I'm used to the °C, In case of K delete the subtract part
 // Add scaled bands to the original image
 return image.addBands(opticalBands, null, true)
 .addBands(thermalBands, null, true);
}

// Function to Mask Clouds and Cloud Shadows in Landsat 8 Imagery
function cloudMask(image) {
  // Define cloud shadow and cloud bitmasks
  var cloudShadowBitmask = (1 << 4);
  var cloudBitmask = (1 << 3);
  var cirrus = (1 << 2);             // I add also mask for cirrus. Suprisely use of Bit 6 (Clear of any clouds) is not working - anybody know why? 
  // Select the Quality Assessment (QA) band for pixel quality information
  var qa = image.select('QA_PIXEL');
  // Create a binary mask to identify clear conditions (both cloud and cloud shadow bits set to 0)
  var mask = qa.bitwiseAnd(cloudShadowBitmask).eq(0)
                .and(qa.bitwiseAnd(cloudBitmask).eq(0))
                .and(qa.bitwiseAnd(cirrus).eq(0));
  // Update the original image, masking out cloud and cloud shadow-affected pixels
  return image.updateMask(mask);
}

// Import and preprocess Landsat 8 imagery
var image = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filterBounds(aoi)
              .filterDate(startDate, endDate)
              .map(applyScaleFactors)
              .map(cloudMask)
              .median()
              .clip(aoi);

// Define visualization parameters for True Color imagery (bands 4, 3, and 2)
var visualization = {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0.0,
  max: 0.2,
};

// Add the processed image to the map with the specified visualization
Map.addLayer(image, visualization, 'True Color 432');
print(image);
// Add the LST Layer to the Map with Custom Visualization Parameters
Map.addLayer(image.select('ST_B10'), {
  min: 19, // Minimum LST value
  max: 43, // Maximum LST value
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ]}, 'Land Surface Temperature');




// // Calculate Normalized Difference Vegetation Index (NDVI)
// var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

// // Define NDVI Visualization Parameters
// var ndviPalette = {
// min: -1,
// max: 1,
// palette: ['blue', 'white', 'green']
// };

// Map.addLayer(ndvi, ndviPalette, 'NDVI Karu')


// Adopted from: https://github.com/panmoravec/GEE/blob/main/LST%20Landsat8


Export.image.toDrive({
  image: image.select('ST_B10'),
  description: 'LST_2023',
  scale: 30,
  region: aoi,
  fileFormat: 'GeoTIFF',
  maxPixels: 3100100100,  // set bigger in case of need
});


// Another version for.... Land surface temperature (LST) Monitoring using Landsat imagery using Google Earth Engine

1) https://code.earthengine.google.com/525ebe138c6d2198a0fd464f7cfc967a 

var LLST = require('users/sofiaermida/landsat_smw_lst:modules/Landsat_LST.js')
var Lan8 = LLST.collection('L8', '2021-1-01', '2021-12-31', geometry, true).select('LST').map(function(image){
  return image.subtract(273.15).clip(geometry).set('system:time_start',image.get('system:time_start'))
})
print("chart");
// print(Lan8)
print(
  ui.Chart.image.series({
    imageCollection: Lan8,
    region: geometry,
    reducer: ee.Reducer.mean(),
    scale: 30,
    xProperty: 'system:time_start'
  }));


Map.centerObject(geometry,10)

Map.addLayer(geometry)

Map.addLayer(Lan8.mean().clip(geometry), imageVisParam4 ,'Dhaka City LST')


====================================================================================


2) https://code.earthengine.google.com/3fdeff0b82b6fc98022a11facd30d269 


// Land surface temperature
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
             .filterBounds(ROI)
             .filterDate('2020-01-01','2022-12-31')
             //.filterMetadata('CLOUD_COVER', 'less_than', 1)
             .median()
             .clip(ROI)
print(landsat8)
Map.addLayer(landsat8, imageVisParam3, 'Satelliteimage with clouds')
Map.centerObject(landsat8,8)


// CLOUD MASK
//Bits 3 and 5 are cloud shadow and cloud, respectively
var  maskL8sr = function(col){
  var cloudShadowBitMask = (1 << 3)
  var cloudsBitMask = (1 << 5)
  
  // Get the piixel QA band
  var qa  = col.select('pixel_qa')
  
   var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
            .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return col.updateMask(mask)
}

// visualization params
var vizParams = {
bands: ['B5', 'B6', 'B4'],
min: 0,
max: 4000,
gamma: [1, 0.9, 1.1]
};

var vizParams2 = {
bands: ['B4', 'B3', 'B2'],
min: 0,
max: 3000,
gamma: 1.4,
};

//load the collection
{
  var col = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
     .map(maskL8sr)
     .filterBounds(ROI)
     .filterDate('2020-01-01','2022-12-31')
     .mean()
     .clip(ROI)
}
print(col)
Map.addLayer(col, imageVisParam2, 'Cloud_Free')

// NDVI CALCULATION

var NDVI = col.normalizedDifference(['B5','B4'])
Map.addLayer(NDVI, {min: -1, max: 0.7, palette: ['blue','white','yellow','red']}, 'Created NDVI')
Map.centerObject(NDVI, 8)

print(NDVI)

// selecting Band10 which is the surface temperature
// thermal band 10
var thermalband  = col.select('B10').multiply(0.1)
Map.addLayer(thermalband, {min:294.388, max:300.514, palette:['yellow','green','blue']}, 'Thermal band')

//finding minimum and maximum NDVI

var min = ee.Number(NDVI.reduceRegion({   //reduce region is used to calculate a statistcal function for each region
reducer: ee.Reducer.min(),
geometry: ROI,
scale: 30,
maxPixels: 1e9
}).values().get(0));
print(min, 'min');
  
var max = ee.Number(NDVI.reduceRegion({
reducer: ee.Reducer.max(),
geometry: ROI,
scale: 30,
maxPixels: 1e9
}).values().get(0));
print(max, 'max')
  
//fractional vegetation
var fv = (NDVI.subtract(min).divide(max.subtract(min))).pow(ee.Number(2))
print(fv,'fv');
Map.addLayer(fv, {},'Fractional Vegetation');

//Emissivity
// you will have to find the emissivity of natural and town
 var a = ee.Number(0.004)
 var b = ee.Number(0.986)
 var EM  = fv.multiply(a).add(b).rename('EMM')
 Map.addLayer(EM, imageVisParam, 'Emissivity')
 
 //FINDING LAND SURFACE TEMPERATURE
 //LST in celsius degree bring -273.15
 //NB: In kelvin don't bring -273.15
 var LST = thermalband.expression(
   '(Tb/(1 + (0.00115*(Tb/1.438))*log(Ep)))-273.15',{
     'Tb':thermalband.select('B10'),
     'Ep':EM.select('EMM')
     
   }).rename('Land surface Temperature');
   var viz =  {min: 22.298638087799763, max:29.22115281845691, palette: [
'040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
'0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
'3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
'ff0000', 'de0101', 'c21301', 'a71001', '911003'
 ]}
   Map.addLayer(LST, viz, 'LST')
   
   
//TIMESSERIES


//print(image_drycollection)

//print(ui.Chart.image.series(LST, ROI, ee.Reducer.mean(), 1000, 'system:time_start'))

// Create the title label.
var title = ui.Label('Shama LandSurface Temperature Map Year of 2022 ');
title.style().set('position', 'top-center');
Map.add(title);

var title = ui.Label('Prepared by Kangah Desmond');
title.style().set('position', 'bottom-center');
Map.add(title);

// set position of panel
var legend = ui.Panel({
style: {
position: 'bottom-left',
padding: '8px 15px'
}
});

// Create the title label.
//var title = ui.Label('Dhaka City LandSurface Temperature Map Year of 2020 ');
//title.style().set('position', 'top-center');
//Map.add(title);
 
// Create legend title
var legendTitle = ui.Label({
value: 'LandSurface Temperature (degC)',
style: {
fontWeight: 'bold',
fontSize: '18px',
margin: '0 0 4px 0',
padding: '0'
}
});
 
// Add the title to the panel
legend.add(legendTitle);
 
// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max-viz.min)/100.0).add(viz.min);
var legendImage = gradient.visualize(viz);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['max'])
],
});
 
legend.add(panel);
 
// create thumbnail from the image
var thumbnail = ui.Thumbnail({
image: legendImage,
params: {bbox:'0,0,10,100', dimensions:'10x200'},
style: {padding: '1px', position: 'bottom-center'}
});
 
// add the thumbnail to the legend
legend.add(thumbnail);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['min'])
],
});
 
legend.add(panel);
 
Map.add(legend);

// Create the title label.
//var title = ui.Label('Prepared By: Mijanur Rahman');
//title.style().set('position', 'bottom-right');
//Map.add(title);

Also Mirza Waleed has a great publication on LST and YouTube explanation on same publication that you should checkout.


Thank you for reading.

Wednesday, January 1, 2025

How to load shapefile into Google Earth Engine (GEE)

 To load a shapefile into Google Earth Engine (GEE), you can use the Code Editor, the Asset Manager, or the command line interface (CLI):

In this guide, I will use the 'Asset Manager' tab and access the shapefile from the Code Editor. Let get started...

Step 1: Go to code.earthengine.google.com and sign in. Click the “Assets” tab at the top left of the screen, and then click on the red “New” button. Select the “Shape files” option from the dropdown menu.


Step 2: In the pop-up window, click the red “Select” button and navigate to where your shapefile is stored. Select all of the files that make up your shapefile (.shp, .shx, .prj, .dbf, .cpg). You can hold the control key to select multiple files at once. Then click “Open” to add the selected files to the upload list. Update the Asset Name in the Asset ID section if needed. Then, click the “Upload” button at the bottom.



Step 3: Navigate to the 'Task' tab, there you see the submitted task in action, you don't have to do anything. Just watch the task will turn blue when the upload is complete.


Step 4: Click on the name of the name of the newly uploaded shapefile under the “Assets” tab to open its details windows where you will all meta data related to it including its "Table ID". The "Table ID" is what we will copy for use in the javascript code editor.



That is it!

Wednesday, December 25, 2024

Preparing Flight Plan for Drone Survey using MissionPlanner

 To conduct a 'Drone Survey' of a site/location, it is required that a proper flight plan should be prepared for the success of the survey mission. Some drones do come with their own software for making the flight plan, however there exist universal tools like the 'Mission Planner' that supports may types of camera on drones thanks to Michael Oborne of Ardupilot.


The 'Mission Planner' software is not only useful for drone surveying, it has many other uses including GNSS RTK (Real-Time Kinematic, Global Navigation Satellite System) processing among others.


Creating the Flight Plan

You will work under the 'PLAN' tab to prepare your 'Flight Plan' for drone surveying.


Switch to the "Plan" tab, then zoom to your study area and position the home point at a suitable location.


Right click and draw polygon over the study area. You may also use the polygon icon to the left for the drawing.



Next, right click and select 'Auto WP' and 'Survey (Grid)'. This will launch a new window where you can select various setting for the drone, camera, forward and side overlaps, images,  etc.

Check the 'Advanced Options' for more configuration settings. When you are satisfied with your settings, click on the 'Accept' button to apply the settings.

Based on the selected study area, type of camera, flying altitude of 120m, camera angle of 360° and flying speed of 20m/s the flight mission has the following statistics:-



Happy Flying!

Sunday, December 22, 2024

Different ways for running python script in QGIS

QGIS was written in  C++, Python and Qt framework.

In this post, I will share with you the common places where you can use or run python script within QGIS. If you already know python and some basics of QGIS, then this will be great exposure to know how to automate different parts of the open source GIS software.

Python in QGIS can be used to interact and automate GIS data in:-

1) Python Console

Plugins >> Python Console (Ctrl+Alt+P)




2) Label Expression



3) Field Calculator



4) Processing Toolbox



5) Plugins



6) Stand alone QGIS Apps



7) Script Runner Plugin



8) Project Event - Open, Close and Save






Thursday, December 12, 2024

Preparing a Flag map of Africa 24 Teams Set For AFCON 2025 In Morocco

 This task consist of two parts, the GIS and Graphic Design parts. To complete the task, I will use two software namely QGIS for the GIS part of the project and GIMP for the graphic design part of the project.

Data Sources:

  • The Map is the world map from QGIS resource data
  • The Flags were made by Hampus Nilsson

Procedures
Step 1: Using the QGIS software, query the world map for African countries and the 24 teams that qualified for the AFCON 2025 as seen below.



Now we know the countries to be designed with flags respectively.

Step 2: Lets make sure the flags we downloaded are available for the 24 countries. The flag source uses two alphabets to name the flags, this could be challenging but fortunately there is a JSON file in the directory the associalte the two two alphabet code with the full name of the countries.


'BW, 'EG', 'ZR', 'SD', 'ZM', 'NG', 'BJ', 'UG', 'TZ', 'GA', 'CI', 'DZ', 'SN', 'CM', 'GQ', 'BF', 'MZ', 'ZA', 'MA', 'TN', 'AO', 'KM', 'ML', 'ZW'

The 'WB_A2' attribute column matches the names of the flags, so that is what we will use the get the names of the country for each flag.

Saturday, December 7, 2024

How to make country flag map in QGIS

Using graphics software like GIMP/Photoshop, Inkscape/Illustrator, we could fill polygons shapes with an image easily, now how do we do the same in QGIS?


Yes, QGIS is not a graphic software yet it has the capability to fill polygons shapes with an image. Let make use of a flag image to fill a country's polygon shape. I will use Nigerian flag and its boundary polygon shape for this demonstration however, you can prety much use any image with any polygon shape of your chosen.

In other words, all we are doing is to insert image in shape in QGIS. Lets see how we can achieve it.

First get the image flag and the vector polygon of Nigeria administrative boundaries.

Now, from the polygon's property window, open symbology tab and select 'Raster Image Fill' as seen below. You can now set path to the image, set the size of the image, set its opacity, etc as you require fit.


One important setting that determine if the image fills the entire polygon layer or each polygon shape in the layer is the "Coord mode" which can either be det to 'object' or 'viewport'.


With some tweaks and tricks, we can achieve what the graphics tool can do in QGIS. Though this kind of map designs are effectively done in the graphics software.

Happy mapping.