texts: Where are the "band names" of an ENVI file once it is imported in R?

texts:
I have a .dat file (ENVI format) with specific band names, as you can see in the snapshot.

I've imported this file in R, but I cannot find the band names in the R object.
The command

names(stack_NDVI)

gives this output

[1] "old_min_identification_VI_smooth_2010.1"  "old_min_identification_VI_smooth_2010.2"
 [3] "old_min_identification_VI_smooth_2010.3"  "old_min_identification_VI_smooth_2010.4"
 [5] "old_min_identification_VI_smooth_2010.5"  "old_min_identification_VI_smooth_2010.6"

that is not correspondent to the band names of the .dat file as I can see in ENVI.

Any ideas to get this information fro the R object?

I've imported the .dat as a RasterStack object.


texts: Create line feature from aspect value in ArcGIS

texts:
I'm investigating the danger of rockfalls into lakes in mountainous areas by using the hydrological tools of ArcGIS 10.7. My problem is that rockfalls don't exactly behave as surface runoff would.

I will try to explain using my crude picture below: My lake is to the right with the light blue outline. The colored areas are slopes with different hazard classifications that possibly drain into the lake.
The dashed line approximately divides these slopes into those, from where a rockfall would hit the lake and from where only runoff water would reach it.

I'm trying to figure out an automated approach to select only the slopes to the right. I have an idea but I don't know if this is possible with ArcGIS:

I would calculate the aspect for every slope pixel and create a line feature according to the aspect value. Then, I would buffer this line with for example 50m and see if it intersects with the lake. This way, I would identify the slopes where more than a certain percentage of pixels somewhat directly point towards the lake.

Now my question:
Is there a way to create a line in arcpy with only a pixel location as origin and an aspect value for direction?


texts: How to use variable field name in arcpy.SearchCursor?

texts:
I have code that says:

fld = fieldName
rows = arcpy.SearchCursor(fc)
for row in rows:
    x = row.fld

How do I use a variable for my field name in Python ArcGIS 10 on Windows?


texts: PostGIS Query: All rows within 1 hour of each other which intersect


texts:
I install Qgis 2.12 but after installation I try to open Qgis and I resive this  three errors.
"The program Can't start because qgis_xxx.dll is missing from your computer. Try reinstalling the programm to fix this problem"
I reinstall Qgis but nothing cahange and the same errors.


texts: Densifying polyline to generate points at set distances using ArcMap?

texts:
I have a polyline, that I wish to create points along using ArcMap 10.2.2.  I used the methodology suggested in Creating series of points along polyline?, but I do not seem to be able to generate the points at the set intervals I desire.  For instance, densifying the polyline using 500m and 250m, and then converting the vertices to points  using the Feature Vertices to Points tool), seems to generate the same number of points (540) for either densifying distance when using the "Both_Ends" option in the Feature Vertices to Points tool:

250m - Both Ends:

500m - Both Ends:

With the same number of points in each generated points layer. And using the option 'All' in the Feature Vertices to Points tool generates way too many points (5250):

My end goal is to have points spaces at every 250m or 500m along each branch of this polyline, finally settling on the polyline to point conversion that generates a reasonable number of values.


texts: Splitting large PostGIS/GeoServer geometry for WMS optmization

texts:
I have a large geometry (MultiPolygon) for a (fictional) world map. It basically spans a large portion of the entire world: it's the difference between a buffered extension of the union of the entire world geometry and that union itself. In other words, a coastal buffer.

The geometry is in a PostGIS table (as one record). I serve this using GeoServer over WMS. Does it make sense to (automatically) split this giant geometry up into smaller ones, or is that completely unnecessary because PostGIS and/or GeoServer use the appropriate (sufficiently fast) intersection queries?

I'm not having a performance problem with this (yet) but the question is nagging at me.


texts: Change domain URL for printURL and createURL

texts:
I have a GeoServer Service enable on my Windows 10 and I install the extension printing.

This URL http://localhost:8010/geoserver/pdf/info.json give me this configuration

{[...],"printURL":"http://localhost:8010/geoserver/pdf/print.pdf","createURL":"http://localhost:8010/geoserver/pdf/create.json"}

I don't know how to change/configure domain name.

I want to change http://localhost:8010/ by https://mydomain.org/


texts: How to display contour labels within QGIS using Maptiler

texts:
How can I display contour labels within QGIS using the Maptiler plugin?
The contour elevation (height) displays 440m in the Identify Result pannel. Can the Layer Properties be reconfigured to display the "height" as a label on the map?


texts: Digital Terrain Analysis - Roughness Index

texts:
Currently I am working on a project dealing with terrain analysis based on a DEM. I came across the term "roughness index" using slope values from that specific DEM. On YouTube I found a video, which showed how to do it in ArcGIS using Focal Statistics:

Creating a slope map based on the DEM.
Using focal statistics to calculate mean, min and max within a 3x3 cell raster for the whole DEM.
Using raster caclulator with this expression:

(mean_elevation - min_elevation) / (max_elevation - min_elevation)

The outcome is called Surface roughness index.

I don't understand why this formula is used.

Why do I have to subtract the minimum elevation from the mean and divide this by the minimum elevation subtracted from the maximum?

Has anyone another (scientific) reference for working with roughness values as standard deviations of slopes?

I have access to ArcGIS 10.4 and QGIS.


texts: I am a beginner in Safe FME and i need to know how to Convert DWG files to JPG Files. Is it possible?


texts: Why is gdalwarp changing value in nodata pixels?

texts:
I'm using gdalwarp to convert MODIS NDVI mosaics from Sinusoidal to epsg:4326. To save space I'm using VRT files until I reach the final product. Here is my workflow.

Mosaic NDVI Modis tiles using pyModis and output as VRT. NoData value (-3000) is kept
gdalwarp the VRT mosaic to epsg:4326, tiff output. Pixels with NoData recieve value 0

My gdalwarp command is:

gdalwarp -t_srs epsg:4326 -tr 0.002 0.002 modis_mosaic.vrt reprojetc_mosaic.tif

Gdalinfo on the reprojected file shows that NoData values are kept (-3000) but nodata pixels in the image are converted to 0.
I'm using gdal 1.11.1 in Ubuntu linux


texts:
Anybody knows how do I clear the tasks list history in Google Earth Engine Editor?
I found some solutions on how to delete tasks with client, but how do I clear the list of tasks in the editor?
I mean this:


texts: When digitizing, objects are not saved well

texts:
I've some problems with digitizing.

If I draw a line or a polygon, the line or polygon is saved in the attribute list, but not visible and its location is not saved either (when I zoom to the object it goes to coordinates 0,0).

Both the project CRS and layer CRS are the same.

Sometimes it happens, sometimes not, and I don't have any idea what I might be doing differently.
If it happens, most of the time QGIS crashes when I try to save the changes.

Thanks,
Silke

I'm using QGIS, Lisboa version


texts: Shaded KML overlays on OSM tiles in QGIS?

texts:
This seems simple - but the solution is escaping me.

I have a series of dynamically generated KML files describing simple, shaded polygons, with each polygon having a slightly different (externally calculated) color and transparency. After import, all features not contained within a polygon would be clipped off.

I'm trying to layer these onto a detailed street map, and thought QGIS might be able to do this.

QGIS doesn't seem to preserve the style attributes from imported KML files. I've tried multiple formats, including running it through an import -> export in Google Earth Pro.

Also, it doesn't appear there is a way to clip the OpenLayers road map tiles outside of the polygon overlay since the OpenLayers isn't showing as a vector layer.

Can QGIS do this or do I need a different piece of software or a web service?


texts:
What is the easiest way to recode data from one geography to another?

For example, if I have population data for one set of polygons, and I'd like to deduce it for different polygons - a bit like a spatial join but taking account of partial overlaps.  For this job I am willing to assume even distribution of population over area.

I imagine you could produce a raster then use zonal statistics to convert back to polygons, but is there a quicker way?


texts: Count rows of each feature class listed in table with ArcPy?

texts:
I have a table in a geodatabase: tableA

A column in that table lists a couple of feature classes that are in the geodatabase.

I want to make a loop on each name, count the number of rows in the attribute table of each feature class and store this value in a column of 'tableA'.

Here is what I've tried:

arcpy.env.workspace =r'C:\Users\user\Desktop\geodatabase\geoA.gdb'
dataset = arcpy.ListTables("tableA") #the table needed is 'tableA'
arcpy.AddField_management("tableA", "rows_count", "LONG") #create field with the count_rows
for i in dataset:
    with arcpy.da.UpdateCursor("tableA",["file_name","rows_count"]) as cur: #the fields that needed
        for row in cur:
            try:
                nrow = arcpy.GetCount_management (row[0]) #count each layer's rows
                row[1].setValue('rows_count', nrow) #fill the column with the values found
                cur.updateRow(row) #apply changes
            except arcpy.ExecuteError as e:
                print(e)

  AttributeError: 'NoneType' object has no attribute 'setValue'

In the result, the rows_count is null for all rows at the creation of the field.

What can be done here?

How the table looks:

file_name  rows_count
file1       null
file2       null
file3 ...   null

Question:

In the main table, how to count the number of rows for each of these files and add store this vlue in a column ?

UPDATE

Wanted output

tableA

file_name  rows_count
file1         0
file2         12
file3         1
...

Where file1,file2, etc are actual feture classes in the gdb and rows_count are the number of rows in each file.


texts: regression model with randomForest in raster package

texts:
I'm learning randomForest classification method in R and I performed land-use classification with randomForest and raster packages. But when I try uesd a regression model produced with randomForest to predict in "raster" with the same operation, there is only a NA value RasterLayer left.

what is wrong, the way I deal with datas or it can't do this in "raster" package


texts: Issue on Updating Geometric Network Weights - Error: -2147210918

texts:
Using ArcGIS 10.4.1, I am trying to add New Weight (Bitgate Type) to an existing Geometric Network but after couple of minutes I am getting this message

  The network weights were not updated
  FDO error: -2147210918
  Cannot register as versioned on this database.

Can you please let me know why this is happening and how I can fix it?


texts:
I'm trying to accomplish something similar to what you have in Google Earth, if you print the map.
You have an area, where you can see what you will get on the paper when you print, which have the right aspect ratio according to the DIN-Format.

I checked the OpenLayers example how to print from here https://openlayers.org/en/latest/examples/export-pdf.html.
But you actually don't get what you see, depend on the DPI I choose I can see map sections beyond the viewport.
What should I keep in mind if I want to implement this feature?


texts: Lining up .asc XY points with DXF using QGIS 2.0

texts: I've exported XY points (nodes of a polyline with attributes) (.asc) from Star APIC and attempted to import it into QGIS as a delimited text file. I cannot however get it to line up correctly with the corresponding DXF file I exported from the same database. My DXF lines up with aerial photographs, so I am sure it is correctly projected. The CSV should be in the same projection (since it comes from the same database). I am considering using the Affine Transformation plugin to rescale and transform. Is this the best solution?


texts: QGIS Python Crashing When trying to Use Layer Functions

texts:
I'm running QGIS 3.10 and trying to build a simple export script for an already configured map file, but I keep running into an issue. My code is below.

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
#from PyQt5.QtCore import *
#from PyQt5.QtGui import *
import os
os.environ['PROJ_LIB'] = "D:/OSGeo4W64/share/proj"
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = "D:/OSGeo4W64/bin/apps/Qt5/plugins"
os.environ['QT_PLUGIN_PATH'] = "D:/OSGeo4W64/apps/Qt5/plugins"
QgsApplication.setPrefixPath("D:/OSGeo4W64/", True)
qgs = QgsApplication([], False)
qgs.initQgis()
project = QgsProject.instance()
print(project.fileName())
project.read("D:\MapWorkingFolder\RanchMaps.qgs")
print(project.fileName())
base_path = os.path.join(QgsProject.instance().homePath())
pdf_path = os.path.join(base_path, "Ranch_Maps_Atlas.pdf")
manager = QgsProject.instance().layoutManager()
layout_1 = manager.layoutByName('ranch_maps_new')
exporter = QgsLayoutExporter(layout_1)
#exporter.exportToPdf(layout_1.atlas(),pdf_path,QgsLayoutExporter.PdfExportSettings())
exporter.exportToPdfs(layout_1.atlas(),pdf_path,QgsLayoutExporter.PdfExportSettings())
#iface.actionExit().trigger()
qgs.exitQgis()

It all "works" aka doesn't throw any errors that I notice, but the exporter command throws a 6 Iterator Error, and if I manually work in Python3 and try to run any function that accesses the layers Python crashes. I've tested this with a more simple map on the same file, and it works.

Example function the crashes this. project.mapLayers()

The map only has 3 layers. 2 PostGIS database layers, and 1 XYZ Tile Layer.

I've searched and tried quite a few things, but without success.


texts: Creating a clean TIN interpolation without "Bleeding" in QGIS

texts:
If you assume I know little to nothing about QGIS 3.10, but have been attempting with some success to create fictional world maps with it, and that I have very limited (OpenScript and Excel) programming skill, perhaps the imprecise technical nature of my question will make more sense!

I am creating Raster files to serve as DEM files.

My process is to create a a set of polygons (within a multi-polygon shapefile layer) for elevation contours.

I have tried starting the coastline contour as elevation 0, elevation 1 and elevation 10

I have used elevation (contour) intervals of 100

I then use Vector/Geometry Tool/ Extract Vertices and get a series of points along each polygon perimeter.

I use the Processing Toolbox / Vector Geometry/ Add Geometry Attributes tool to create x,y,z coordinates for each point

I have then gone into the attribute table, selected all of the points, used the multi edit tool and set zcoord = to elevation

I then have either created contours from the data, or just used the points themselves with similar results, and completed a TIN interpolation.

The raster files are sufficient to create a decent 3-D map, but the triangulation plays havoc with my coastline (the 0,1 or 10 elevation contour) as shown in the photo.

You can see the triangles departing the coastline, which produces a raster that doesn't really show the coast at all.   I have added a second photo to demonstrate and you can see the triangles jutting off the coastline in the raster file

I am certainly committing some error in this process and wondering if anyone can direct me to a tutorial or provide some suggestions for how to create a more precise coastline.

I have not yet attempted creating a -100 contour around the island and then a 0 elevation coastline.

Here it is after implementing Gabriel De Luca's suggestion about adding points where sections of the polygons were excessively long unbroken lines.


texts:
I am trying to create a buffer around several points using the function circles in the R package dismo. The problem arises with points that are close to the 180 or -180 longitude. In these cases the buffer is not a circle but a band that connect both sides of the map extension. Probably I am missing some basic option, but I am wondering if there is any way to avoid this issue. Below an example code and its result.

library(dismo)
library(dismo)
xx<-cbind(c(-175.20, -106.65,-103.97,-17.76),c(-21.13,35.08,36.78,28.65))
pp<-as.data.frame(rep(1,4))
zz<-SpatialPointsDataFrame(xx,pp)
cc<-circles(zz,1000000,lonlat=TRUE)
plot(cc)


texts: QGIS 2.14.2 Rotating map quick map service

texts:
I am currently using an Esri basemap.
Standard built in functionality of QGIS 2.14.2 allows for the rotation of vectors, rasters etc. When I use rotate, all the vectors I added personally rotate just fine.
However, rotation completely screws up the view of my basemap (esri topo. however the same holds for any other map I loaded in via the "Quick map service"). Whenever I rotate, QGIS does something to the basemap, tranforming the entire look of the map.

How can I rotate my basemap with it changing completely?

p.s. even if I try in map composer, the same thing happens.

----------UPDATE----------

Here are some pictures to help clarify to issue...

The funny thing is, that the shape files, do rotate correctly.
But the underlying basemap, gets completely warped. How to prevent this


texts: Find point on line with given nearest point

texts: I have a line and point layer. Points are not exactly intersect with line layer. With respect to existing point, how to find the new points which are exactly on the line ? consider buffer distance is provided.


texts: Synergis Weboffice vs Geocortex Essentials?

texts:
What are the main differences between Synergis Weboffice and Geocortex Essentials?

I know they are both used for creating web-based mapping application for end users and they depend on ArcGIS server for working but are there any differences between the tools and features they provide?


texts: Create 2nd UI Window for Qgis 3.0 Plugin

texts:
I have built a Qgis 3.0 Plugin which was created by QGIS Plugin Builder. The Plugin has two separate options which work perfectly (see Figure).

However, I want to attach two completely separate UIs for the two options. For the first one I already used the .ui file from the Skeleton, created by the Plugin Builder.

My Question is now: How can I create a second UI Window for the second option. It should only show two QFileDialog Items where the user can select two import directories and a run-Button.

I assume I must edit the parent attribute in the initGui method and set up another .ui instance?

For completeness a code snippet from my initGui method:

def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    export_icon = ':/plugins/export_rf3d/export.png'
    self.add_action(
        export_icon,
        text=self.tr(u'Export vector files to rockyfor3d ascii files'),
        callback=self.run,
        parent=self.iface.mainWindow())

    """Import rf3d generated files and style them"""
    import_icon = ':/plugins/export_rf3d/import.png'
    self.add_action(
        import_icon,
        text=self.tr(u'Import Files from rockyfor3D output directory'),
        callback=self.runImportRockyFor3DFiles,
        parent=self.iface.mainWindow())


texts:
I'm trying to obtain Sentinel-1 images in Google eart hengine and to have 3 different imagecollections:
1. HV
2.HH
3.VV
 I have written the following code but whenever I check my image collection it has 0 elemnts. I know that there are images on the chosed period. what can  be the reason? where is my mistake?
 I'ts important to mention that in the following code, everytime I used "print" in order to check my image collectio,  I have gotten 0.

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter by metadata properties.
var dataset = sentinel1
  // Filter to get images with VV and VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'HH'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'HV'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .filterBounds(geometry);

// Filter to get images from different look angles.
var descending = dataset.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

print(descending);
//Filter dates
var dateDesc = descending.filterDate(ee.Date('2018-10-01'),ee.Date('2019-10-08'));

print(dateDesc);

//clip the dataset according to the geometry
var clippedCol=dateDesc.map(function(im){
  return im.clip(geometry);
});

var HV=clippedCol.select('HV');
var HH=clippedCol.select('HH');
var VV=clippedCol.select('VV');
print(HV);

I have also tried this code but again , got 0 elements:

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter VH, IW
var HV = sentinel1
  // Filter to get images with VV and VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'HV'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  // reduce to VH polarization
  .select('HV')
  // filter 10m resolution
  .filter(ee.Filter.eq('resolution_meters', 10))
  .filterBounds(geometry);
// Filter to orbitdirection Descending
var HVdescending = HV.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
// Filter time
var HVdate = HVdescending.filterDate(ee.Date('2015-08-01'), ee.Date('2019-10-08'));

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter HH, IW
var HH = sentinel1
  // Filter to get images with VV and VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'HH'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  // reduce to VH polarization
  .select('HH')
  // filter 10m resolution
  .filter(ee.Filter.eq('resolution_meters', 10))
  .filterBounds(geometry);
// Filter to orbitdirection Descending
var HHdescending = HH.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
// Filter time
var HHdate = HHdescending.filterDate(ee.Date('2015-08-01'), ee.Date('2019-10-08'))

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter VH, IW
var VV = sentinel1
  // Filter to get images with VV and VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  // reduce to VH polarization
  .select('VV')
  // filter 10m resolution
  .filter(ee.Filter.eq('resolution_meters', 10))
  .filterBounds(geometry);
// Filter to orbitdirection Descending
var VVdescending = HV.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
// Filter time
var VVdate = VVdescending.filterDate(ee.Date('2015-08-01'), ee.Date('2019-10-08'))

print(VVdate);
print(HVdate);
print(HHdate);


texts: Renaming legend items in ArcGIS for Desktop?

texts:
Can someone please tell if it is possible to rename legend items in the map layout in ArcGIS for Desktop?

If yes, how can I do it?


texts: PyQt5: QFileDialog triggered as many times as plugin is opened

texts:
I am new to QGIS and I'm building a plugin on it. The problem I'm facing, I believe, concerns mostly PyQt5.

So, in my plugin I'm trying to open a file through getOpenFileName() of QFileDialog. I noticed that the dialog for opening the file is triggered as many times I trigger the plugin without reloading (using Plugin Reloader). Putting it simply, the dialog for opening the file pops up as many times in a row as I test the plugin without reloading it. The first time user has to input file name once, the second time twice, third time thrice and so on. I want this dialog to pop up just once and get the input from the user.

I'm using Plugin Builder. Here's my code for reference:

def loadData(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop

    self.dlg.lineEdit.clear()
    self.dlg.toolButton.clicked.connect(self.select_output_file)

    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # ...

def select_output_file(self):
    filename = QFileDialog.getOpenFileName(self.dlg, "Select file", "")[0]
    self.dlg.lineEdit.setText(self.tr(filename))

Here, lineEdit is the name for text field where the path of the file being opened appears, toolButton is the button that is clicked to trigger this action and dlg is the instance for main dialog box of the plugin.


texts: postgis - managing single dataset, multiple coordinate systems

texts:
I am setting up a forest inventory database that aims to store information about forest plot data from various countries. The current project is in Mozambique, others are in Laos, etc.

Each project provides their forest data in different coordinate systems. Other than that the plot data (fields) is identical (nearly), so it makes sense to store it all in one table.

I have come across the problem today where I want to view multiple project data in Mozambique. The CRS from one project is supplied in WGS84, the other UTM 36S. The geometry table has set the CRS for plots to WGS84, so when I load in QGIS I can't see the plots in UTM.

What is best practice?

To choose a single CRS (WGS84) for all plot data?
To store plot data from different countries in different tables?
Or depending on which record I am viewing QGIS can dynamically use the relevant projection?


texts: Capturing Encrypted String from ESRI Python Toolbox Parameter

texts:
I created an ESRI Python toolbox and setup one of the parameters of the toolbox as an encrypted string as the parameter is a location for the end user to enter a password.

The problem I am encountering is getting that string un-encrypted from the parameter using the only method I have ever used parameters[0].valueAsText and the only method I could find in the ESRI documentation. This method keeps the string encrypted. Subsequently, I encounter errors when trying to use this parameter in the rest of the script. I'm using ArcGIS Desktop 10.3

 def getParameterInfo(self):
    """Define parameter definitions"""

    param0 = arcpy.Parameter(displayName = 'AGOL Username',
    name = 'AGOL Username',
    datatype ='GPString',
    parameterType ="Required",
    direction ="Input")

    param1 = arcpy.Parameter(displayName = 'AGOL Password',
    name = 'AGOL Password',
    datatype ='GPEncryptedString',
    parameterType ="Required",
    direction ="Input")

    param2 = arcpy.Parameter(displayName = 'AGOL Group ID',
    name = 'AGOL Group ID found in the URL, when viewing group in AGOL',
    datatype ='GPString',
    parameterType ="Required",
    direction ="Input")

    param3 = arcpy.Parameter(displayName = 'Output Folder',
    name = 'OutPut Folder to Store Zipped Downloaded Feature Services',
    datatype ='DEFolder',
    parameterType ="Required",
    direction ="Input")

    params = [param0, param1, param2, param3]
    return params

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    return

def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

def execute(self, parameters, messages):
    """The source code of the tool."""
    import os #,tempfile, shutil
    #import zipfile
    import json #, re
    #import uuid
    import  time #datetime,
    import urllib,urllib2
    #import sys

    username = parameters[0].valueAsText
    password = parameters[1].valueAsText
    groupid =  parameters[2].valueAsText
    outputLocation = parameters[3].valueAsText


texts: Spatial join geotagged images to preexisting point shapefile?

texts: I have a layer of geotagged images (points with jpg thumbnail attachment), and another layer of different points (shapefile with various attributes). I would like to perform a spatial operation in which the jpg thumbnail from the geotagged image layer gets added to the attribute table of the closest point in the other layer. I've tried spatial join but I'm not sure that I am doing it correctly. Any ideas?


texts: Conversion between semicircles and latitude units

texts:
I've recently come across the semicircle unit that is used to describe a latitude and longitude using the maximum precision allowed by a 32-bit number. The conversion here is given as

degrees = semicircles * ( 180 / 2^31 )
semicircles = degrees * ( 2^31 / 180 )

However it doesn't explain if this should be used differently for latitude since the maximum allowed number of semicircles is 2^31, and if you substitute this into the equation you get the answer "180 degrees", but latitude only goes up to 90.
Would a different conversion between latitude be used instead? Such as:

degrees = semicircles * ( 90 / 2^31 )
semicircles = degrees * ( 2^31 / 90 )

If so then you wouldn't be using the 'semicircle' as a unit for latitude, more like a quarter circle in this case.
Has anybody else come across this issue?


texts:
I'm working with 2.14, plugin QGISCloud (free version) and for a shapefile (point) I want to pop-up a picture of my vacation and share it with others. It works offline with the code in the HTML box / expressions:

[%CONCAT('<img src="', "Loc", '" width="400"></img')%]

"Loc", stands for name and folder direction on my local drive. Backslash like: / (backward-sloping) not \(upward-sloping).
But my question is: How to do it online with QGISCloud? Are there alternatives (with other plugin's).


texts: How edit a table for get a new one in ArcGIS?

texts:
I have a shapefile and y want to transform the attribute table (image 1) in other table (image 2) but I only want to get some fields and values.

From the atribute table I want to transform FID_malla field as rows and the GROUP fields values as a columns in unique values. Then the table values are filled with Shape_Area field values.

I use R language to do it but, there is any tool in ArcGIS to do that? or I have to do it using ArcPY (UpdateCursor)?

(image 1)

(image 2)


texts: how to get the value from ui.Chart.image.series

texts:
I print a chart in the Earth Engine code editor using ui.Chart.series. I would like to extract the values of this chart using the Earth Engine API in Python to produce my own chart with matplotlib.
here l8 is any ee.imageCollection having 'cloud and 'NDMI' band and point is a ee.Geometry.Point.

var buffer = point.buffer(30);
print(ui.Chart.image.series(l8.select(['cloud', 'NDMI']), buffer, ee.Reducer.mean(), 30));

Is there a function to extract the data as a ee.list or a ee.dictionary or can I found somewhere what is under the hood of this function ?


texts: How can I import assets of type folder shared by my co-worker?

texts: I am very new to Google Earth Engine. My coworker shared his asset folder with me but as i open the link it gives me an error that "Unable to import asset of type folder".


texts: Does ArcGIS Directions API take Traffic into consideration when calculating Driving Time?

texts:
I'm trying to calculate the Estimated Time of Arrival (ETA) between 2 locations and I need to take traffic into consideration.

I found this API call but it isn't clear if it includes traffic in the time estimates.

For example, the RouteService response below has "totalDriveTime"

"directions": [
        {
            "routeId": <routeId>,
            "routeName": "<routeName>",
            "summary": {
                "totalLength": <totalLength>,
                "totalTime": <totalTime>,
                "totalDriveTime": <totalDriveTime>,
                "envelope": {
                    <envelope>
                }
            },


texts:
There is no problem in running this code in qgis 2.x but it through error when I tried in qgis 3.

def _calculate_means(self, polyg_layer_files, raster_layer_files, attribute_names):
    # polygon_layer_files: list of absolute path polygon layer files
    # raster_layer_files: list of absolute path raster layer files
    # attribute_names: name prefix of newly calculated attributes (means from raster to polygon)
        for poly in polyg_layer_files:
            for rast, attr_name in  zip(raster_layer_files, attribute_names):
                assert os.path.isfile(poly), "File {} does not exist.".format(poly)
                assert os.path.isfile(rast), "File {} does not exist.".format(rast)
                poly_object = QgsVectorLayer(poly, 'zonepolygons', 'ogr')
                zoneStat = QgsZonalStatistics(poly_object, rast, attr_name, 1, QgsZonalStatistics.Mean)
                check = zoneStat.calculateStatistics(None)
                assert check == 0, "zoneStat.calculateStatistics(None) returned non-zero value... check input layers!"

I am not able to figure what is the actual problem here. The error that I got is:

File "/home/work_dir/my_code.py", line 98, in _calculate_means
    zoneStat = QgsZonalStatistics(poly_object, rast, attr_name, 1, QgsZonalStatistics.Mean)
TypeError: QgsZonalStatistics(): arguments did not match any overloaded call:
  overload 1: argument 2 has unexpected type 'str'
  overload 2: argument 1 has unexpected type 'QgsVectorLayer'


texts: Merging entities to existing feature class?

texts: Do you have any idea how I can merge entities to an existing feature class that already has the same structure and not in an existing feature class as the ArcGIS tool proposes?


texts:
I was wondering which axis X Y and Z point which way for PostGIS BOX3D-Objects.

I thought Z was up, X was right and Y was "back", which matches this drawing: https://gis.stackexchange.com/a/214938/120335

However, looking at the code here, I found a comment calling the side with the maximum Y front, which would mean the Y axis would point the other way.

So, is one way wrong and one way correct? Are both correct in their own ways?


texts: On load event for Leaflet proj4leaflet projected geoJSON layer

texts:
Can't find proper event for vector layer like on "load", to be fired when this layer is fuly rendered. I add vector features to a layer in Leaflet using proj4leaflet from JSON object like this

   var geojsonLayer = L.Proj.geoJson(features, {
        'pointToLayer': function (feature, latlng) {
            return L.marker(latlng);
        }
    }).addTo(map);

And this methid does not work

geojsonLayer.on("load", function () {
        console.log(1);
    });

I need a loadend or renderend event for this layer. Is there any?


texts: Exporting an image for visualization (Earth Engine Python API)

texts:
In a snippet from the documentation I can view elevation data in the code editor with the following:

var dataset = ee.Image('JAXA/ALOS/AW3D30/V2_2');
var elevation = dataset.select('AVE_DSM');
var elevationVis = {
  min: -9999,
  max: 15355,
  palette: ['0000ff', '00ffff', 'ffff00', 'ff0000', 'ffffff'],
};
Map.setCenter(136.85, 37.37, 4);
Map.addLayer(elevation, elevationVis, 'Elevation');

The palette and min/max values are defined in a layer inside a map. How do I access them in Python?

So in the following code for export to drive, how would I make the DEM data be something useful? Right now its just a black mess since I can't control the visualization.

dataset = ee.Image('JAXA/ALOS/AW3D30/V2_2').select('AVE_DSM')

task = ee.batch.Export.image.toDrive(image=dataset,
                                    description= 'Image',
                                    folder='BatchTest',
                                    region=[[-120, 35], [-119, 35], [-119, 34], [-120, 34]],
                                    scale=30,
                                    fileFormat='GeoTIFF')
task.start()
print(task.status())


texts: Measuring distance of image pixel from different neighboring polylines?

texts:
I'm going to measure the distance between a given image pixel and a series of neighboring polylines, consisting of 15 different types of roads. In other words, I'm looking for a solution with which I can measure the distance of the pixel from all the 15 types of roads, not only the nearest one.

How can I do it in ArcGIS?


texts: Does DotSpatial Reproject work for non-global to global projections?

texts:
As a small part of a project i need to convert shapefile geometry from OSGB1936 (EPSG:27700) to WGS1984/Pseudo (EPSG:3857) and output as wkt. I found the DotSpatial library which looked promising and shoehorned this;

public void button1_Click(object sender, EventArgs e)
        {

        var pathin = "G:\\SearchArea_region.shp";
        var pathout = "G:\\wkt.txt";

        using (var fs = FeatureSet.Open(pathin))
        {
            var writer = new WktWriter();
            var numRows = fs.NumRows();
            for (int i = 0; i < numRows; i++)
            {
                var shape = fs.GetShape(i, true);
                var geometry = shape.ToGeometry();

                double[] xyArray = new double[geometry.Coordinates.Count() * 2]; //START Reprojection
                double[] zArray = new double[1];
                int counterX = 0;
                int counterY = 1;
                foreach (var coordinate in geometry.Coordinates)
                {
                    xyArray[counterX] = coordinate.X;
                    xyArray[counterY] = coordinate.Y;
                    counterX = counterX + 2;
                    counterY = counterY + 2;
                }
                var EPSG27700 = new ProjectionInfo();
                var EPSG3857 = new ProjectionInfo();
                EPSG27700 = KnownCoordinateSystems.Geographic.Europe.OSGB1936; //OSBNG
                EPSG3857 = KnownCoordinateSystems.Geographic.World.WGS1984; //Pseudo Merc
                Reproject.ReprojectPoints(xyArray, zArray, EPSG27700, EPSG3857, 0, (xyArray.Length / 2));
                counterX = 0;
                counterY = 1;
                foreach (var coordinate in geometry.Coordinates)
                {
                    coordinate.X = xyArray[counterX];
                    coordinate.Y = xyArray[counterY];

                    counterX = counterX + 2;
                    counterY = counterY + 2;
                }
                geometry.GeometryChanged(); //END

                var wkt = writer.Write((Geometry)geometry);
                MessageBox.Show(wkt);
                System.IO.File.WriteAllText(pathout, wkt);
            }

        }

    }

From a question elsewhere on this site. However, this outputs only "infinity/infinity" coords or does no reprojection at all. I have tried to find DotSpatial tutorials but links back just seem to go back to the archived site landing page with no accessible information.

Either, the code is rubbish but i don't know enough to fix it or the re-project is failing because OSGB isn't a global system.


texts: Where can I find USGS SRTM 3 index grid or what procedure should I employ to select tiles?

texts:
I would like to download all USGS SRTM 3 tiles for Ontario, Canada from the USGS site: http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/

Where can I obtain the index that will assist me in querying the tiles I need to download. Alternatively, is there another easier way to identify the tiles?

When I was downloading the tiles for Manitoba, Canada I found a range of files then just downloaded the individual files from the site. Is there a similar range of files for Ontario?

Thanks


texts: Loading KML WMS Network Link through Google Earth API

texts:
I'm trying to add a KML network link with the Google Earth API using the same function that I have used for regular KMLs and KMZs. Sources suggest to create a KML with a "Ground Overlay" to connect to the WMS which I have done and displays perfectly well within Google Earth. My problem lies when I attempt to add the same KML with the Google Earth API in my application.

My function to load KMLs & KMZs to the map:

var geoNat1 = new Array(11);
var contNat1 = 0;

function nacional(tip, nombre, numer, obj) {
    contNat1 = numer;
    $j('.cargando').show();
    console.log(nombre);
    if (geoNat1[contNat1] == null) {
        geoNat1[contNat1] = ge.createNetworkLink("");
        geoNat1[contNat1].setDescription("NetworkLink se abre para el contenido recuperado");
        geoNat1[contNat1].setName("Abrir NetworkLink");
        geoNat1[contNat1].setFlyToView(true);
        var link = ge.createLink("");
        link.setHref("http://" + servidor + "/" + directorio + "/carto_ref/" + nombre);
        geoNat1[contNat1].setLink(link);
        ge.getFeatures().appendChild(geoNat1[contNat1]);
    } else {
        if (obj.checked == true) {
            geoNat1[contNat1].setVisibility(true);
        } else {
            geoNat1[contNat1].setVisibility(false);
        }
    }
    $j('.cargando').hide();

}

My KML:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<GroundOverlay>
    <name>nowCOAST - Tropical Cyclones</name>
    <drawOrder>1</drawOrder>
    <Icon>
        <href>http://nowcoast.noaa.gov/wms/com.esri.wms.Esrimap/wwa?SERVICE=WMS&amp;REQUEST=GetMap&amp;FORMAT=image%2Fgif&amp;TRANSPARENT=TRUE&amp;STYLES=&amp;VERSION=1.1.1&amp;LAYERS=NHC_TRACK_POLY%2CNHC_TRACK_LIN%2CNHC_TRACK_PT%2CNHC_TRACK_WWLIN%2CNHC_TRACK_PT_72DATE%2CNHC_TRACK_PT_120DATE%2CNHC_TRACK_PT_0NAMEDATE%2CNHC_TRACK_PT_MSLPLABELS%2CNHC_TRACK_PT_72WLBL%2CNHC_TRACK_PT_120WLBL%2CNHC_TRACK_PT_72CAT%2CNHC_TRACK_PT_120CAT&amp;WIDTH=1500&amp;HEIGHT=637&amp;SRS=EPSG:4326&amp;BBOX=-118.867172,14.532850,-86.703392,32.718620</href>
        <refreshMode>onInterval</refreshMode>
        <refreshInterval>3600</refreshInterval>
    </Icon>
    <LatLonBox>
        <north>32.71862</north>
        <south>14.53285</south>
        <east>-86.70339199999999</east>
        <west>-118.867172</west>
    </LatLonBox>
</GroundOverlay>
</kml>


texts:
I would like to execute the Curves And Lines add in from my Python script.

Is there any interface for working with add-in inside of a Python script?


texts: In QGIS is it possible to apply an expression to a join

texts:
I am trying to join by 2 text fields and there is a lot of mismatch due to CaSe and space issues... can it be run so that it's

lower("field1")=lower("field2")

I otherwise have to create a new field with this and as it's a repeated process I'd rather not.

If Virtual Layer is the answer -what's the SQL query to run?

SELECT from table1, table2
WHERE lower(table1.field1)=lower(table2.field2)
GROUP By field1

Ideally I'd like to use RE and remove extra spaces at the end as well.


texts: Converting raster with attribute table to ASCII using ArcGIS Desktop?

texts:
I have a vegetation GNN data downloaded from here https://lemma.forestry.oregonstate.edu/data/structure-maps
It comes as a grid raster file with an attribute table, but I need to convert it to a raster file in ASCII format based of certain columns.

Is there a way to do that in ArcGIS Desktop?

I tried Delete Raster Attribute Table Tool, it runs successfully, but does nothing.


texts: Replace non-ASCII characters in QGIS Field Calculator

texts:
I want to replace all non-standard characters from fields in a dbf table in QGIS.

From https://stackoverflow.com/questions/20078816/replace-non-ascii-characters-with-a-single-space

I see we can use something like

re.sub(r'[^\x00-\x7F]+','_',  "New_Name" )

But I am not sure how to use this in the the Expression or Function Editor module.

So I have a column called "New_Name" which I want the !@#$%^&*()~`+ etc signs removed and replaced with _

Can RE be used or do we need to specify each value to replace using CASE WHEN or Replace -the solution in RE is better if it can be done in QGIS Field Calculator using the regexp_replace expression?


texts: How to create a cloropleth world map in ArcGIS 10.2

texts:
I am new to ArcGIS and I am trying to create a choropleth map to display ratio values for different countries.

I believe that my basemap will need to have a country layer so that the shading of country polygons will be possible but I don't understand how to link this with my data.

Will the world topographic map have an attribute table for the country layer in which I can add a new column containing my ratio values or do I need to create and import an excel file contains this information? If I need to create an excel sheet how do I link this so that the country polygons fill to represent the ratio values?

I have found some useful tutorials to show how I can view layer properties and select quantities from the symbology tab however I don't understand have to get my data to this stage.

Please can someone point me in the right direction

Thanks!


texts: Overpass API Python wrapper (incomplete polygon error)

texts:
I am trying to use the following query to get all buildings from the given area in OSM with the Overpass Api, but I am receiving an "Received corrupt data from Overpass (incomplete polygon)" Error. Where is the error in my query?

import overpass
api = overpass.API()

query = """
[out:json][timeout:25];
(way["building"](47.35436014253597,8.490386538763183,47.39245772996116,8.57950135269201);
relation["building"](47.35436014253597,8.490386538763183,47.39245772996116,8.57950135269201);
);
out body;
>;
out skel qt;
"""

response = api.get(query)


texts: Calculating percentage of cells which are lower than site in Grass GIS?

texts:
In GRASS GIS:

I have a shapefile with 300 hundred of sites (with height values for each point), and want to calculate the percentage or number of cells (of a given DEM) which are lower than the site's value.


texts: DSM void filling results into problems with pixels (using ArcMap)

texts:
Before writing here I exhausted any online resources, blogs, videos etc I could find. Here's my problem, with photos and a thorough explanation.
I have a very high res DSM (5m) that has holes (voids). I need to keep and use this surface and I am searching for a way to make it adequate for the usage I want it, that is for hiking maps. Searching for a replacement of lower res is not a solution for my case. Also, I am currently using ArcMap and do not spare the time to switch to ArcGIS Pro right now.
Trying to fix this, I came up with the following workaround that still produces some visible inadequacies (see pictures 4 and 5):

I smoothed the DSM 5m with focal statistics in order to smooth the edges of the voids
I used a DTM 30m (ALOS) which I clipped at the same extent and then resampled at 10m.
I combined my DSM 5m and the DTM 10m with map algebra (raster calculator) with the following expression: Con(IsNull("DSM_5m"),"DTM_10m","DSM_5m")
I examined the result and even smoothed it with focal statistics (circle, 3px radius), but it is still inadequate.

My problem is: at the edges of the void fill, the merging/stitching is very visible in surfaces like hillshade or slope and derivatives like contour lines, a result which is not acceptable for hiking maps.
In the photos below I pinpointed a few pixel areas where the problem is very visible.
My question is: is there a way to fix this, i.e. smooth the edges before or after combining the two rasters in order to fill the voids?
Perhaps I have to start over with another procedure, in case I am doing something wrong or not really helpful (like combining the two rasters)
The rasters I'm combining. The DTM 30m (only visible here in the holes as grayscale) fills the voids of the DSM 5m

The rasters I'm combining. Hillshade for DSM 5m

The rasters I'm combining. Hillshade for both rasters.

Showing the problem at the edges of the void fills

Showing the problem at the edges of the void fills #2


texts: ArcGIS Web Mapping Service Extract By Mask Give 99999 Error

texts:
The USDA used to provide four-band NAIP imagery as TIFF tiles for each township (6x6 miles) in the USA.  The 2013 data only available as an ArcGIS web-mapping service.  Not being a big fan of streaming data I would like to download these data for archiving and speedy rendering.  If I zoom to a township size extent and try Data-Export Data (limited by data frame extent) I get an error:

Output raster dataset exceeds the size limitation (rows and columns) defined on the server. Export a smaller area by reducing the number of rows and columns...

I can zoom into a section (1 mile x 1 mile) and execute a data export without error.

Being clever I thought I would automate the extraction of section tiles with Python.  As a test I have added a township-range-section feature class and selected a single section.  I then tried running Extract by Mask using the imagery and the selected section.  This returns the following error:

ERROR 999999: Error executing function.
The requested image exceeds the size limit.
Failed to execute (ExtractByMask).

Since I can download a 1 mile by 1 mile tile manually this extent must be OK so  why does the Extract by Mask function fail?  Is there something about the WMS tile that cause the geoprocessor to fail to recognize the selected section when performing the Extract By Mask?


texts: Leaflet Control Search and GeoServer WFS

texts:
Maybe someone could help me how to use plugin Leaflet Control Search with WFS from GeoServer. The control is added to the map. But when I try to search a single feature with this plugin, the WFS response with full list of all features. How should I use the Leaflet Control Search to get ony one, searched feature from WFS. I'm quite new in Leaflet and GeoServer. Please, help. Here is what I made:

mapid.addControl(new L.Control.Search({
   url: "http://localhost:8000/geoserver/przejazdy/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=przejazdy:lokalizacje&outputFormat=application/json"
   propertyLoc: ['wsp_lat', 'wsp_lon'],
   propertyName: 'nr_id',
   zoom: 16,
   markerLocation: true
}));


texts: Joining lines when direction is not known

texts:
I am struggling with this problem. I have a polyline made up of several segments. Each segment is made up of points, each point identified in 3D space by coordinates and elevation. Together, if plotted, the segments make up a more or less continuous line (there could be breaks between segments), but the segments themselves are not consecutive, nor do the points of all segments follow the same direction of travel.
The question is: how can I create, using Python preferably, a single line from these non-consecutive segments so that I can measure the distance between points and the total length of the line.
I don't even know which of the segments is the first one or the last one in the line, but somehow I have to put them in the right sequence and make sure they all point in the same direction so I can measure them.
thank you for any assistance. I'll be happy to supply additional info, data, etc..
I stress I am not asking for the actual python code (not that I would refuse it...), just the strategy.
Bob


texts:
There is a new accurate projection for world maps ► Equal Earth.

I know its proj4 but I can't add it in QGIS (see below). Does it need some kind of QGIS update?


texts: Can't Register table with Geodatabase - ArcGIS 10.5, PostgresSQL 9.5

texts:
I have made connection to a postgresql database in ArcGIS Catalog from Database connections . pgAdmin shows the database owner is 'postgres' so i have made the connection with user 'postgres'.

I want to further register a table in this database to Geodatabase (to add OID field) which i am not able to do. When I right-click the table and go to Manage, everything is greyed out and not available.

When i use GP tool 'Register with Geodatabase' it gives error -
ERROR 001050 Either registered with geodatabase already or cannot open the dataset.

Documents on ArcGIS resource say - 'Be sure you connect as the owner of the table.'

What does it mean here ? Does it refer to the owner of the database in postgresql or .. ?


texts:
What is the method to close labels and features on satellite images of Google maps?

NOTE: I don't want to use classic version of google map or use of history of google map.

What i want to remove or close, are shown in given example image.


texts: Filter sentinel 2 images that have empty areas

texts:
When trying to download a clipped area from the sentinelhub API python library, it returns images from 2 different orbits: One that passes right on top (Right pic), and the image is full, and a orbit next to (Left pic), but its image is cut due to deformations.

Checking the metadata in the sentinelhub website, there is a parameter called no_data_pixel-percentage, but that parameter is not included in the metadata returned by the API. Does anyone have found a solution to filter the good ones?


texts:
I'm pretty new to working with shapesfiles. I'm trying to create a spatial index of addresses in NYC. I completed it using an rtree but have only done it for Manhattan. To complete the other counties in NYC I'd like to combine the shapefiles for each county. The source I'm using for the shapefiles in the NYC Pluto file (http://www.nyc.gov/html/dcp/html/bytes/dwn_pluto_mappluto.shtml#mappluto). It looks like there's isn't a single shapefile for the whole city. Rather there are five shapefiles, one for each county.

Can anyone please point me in the right direction to combine the shapefiles so I can run the rtree code on a singe shapefile, rather than 5.

I'm using Python.


texts: arcpy.da.SearchCursor not work when run the whole scripts in IDLE

texts:
I'm using ArcGIS 10.3.1. The code is running in IDLE of winpython 2.7 64bit.
Here is a very simple code which intends to read a shapefile.

import arcpy
from arcpy import env
env.workplace = "E:/"
shpfile = "landuse.shp"
polygonlist=[]
with arcpy.da.SearchCursor( shpfile,("SHAPE@") ) as cursor:
    for line in cursor:
        polygonlist.append( line[0] )

When I run the script file with IDLE, it failed with "RuntimeError: cannot open 'landuse.shp'. " of the "with arcpy.da.SearchCursor( shpfile,("SHAPE@",col) ) as cursor:"

However, when I type the script lines one by one into the IDLE shell. It worked just fine.

When I run the script file in Spyder. It ended with the same error at the first time. But I just press F5 for another try, it worked fine.

I really want to know why.


texts: How to add labels to shapefile point layer in R?

texts:
I have plotted a shapefile containing points in R, and I would like to add labels like : point 1, point 2 and so on..) to the plot.

Is this possible?


texts:
I am generating a image consisting of 3 bands and would like to generate a thumbnail from it. Unfortunately the thumbnail is not displayed, only a placeholder is shown.

My script is the following (Script in Earth Engine Code Editor)
var europe_bbox = ee.Geometry.Polygon(
    [[[-12.222935630617995, 47.34321448670071],
      [-12.222935630617995, 34.021613745857614],
      [7.288783119382005, 34.021613745857614],
      [7.288783119382005, 47.34321448670071]]], null, false);

// Make a black background image
var background = ee.Image(0).visualize({palette: '000000'});

// Make a grey Europe outline
var addID = function(feature) { return feature.set({ID: 1});};

var europe = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
           .filter(ee.Filter.eq('wld_rgn', 'Europe'))
           .map(addID);

// Create an image from the feature collection
var europeImg = europe.reduceToImage(['ID'], ee.Reducer.first());

//Make the tree loss layer
var treeLoss = ee.Image('UMD/hansen/global_forest_change_2019_v1_7').select('lossyear');

// Blend background - dark europe outline - tree loss layer
treeLoss = treeLoss.visualize({ palette: ['yellow', 'red'], min: 0, max: 19});

var composite = background.blend(europeImg).blend(treeLoss);

print(ui.Thumbnail(composite,
 {format: 'PNG', dimensions: 2000, region: europe_bbox, crs:'EPSG:3035'}));

Does anyone know why this happens and what I may be doing wrong?


texts: QGIS point displacement of only unique point classes/symbols

texts:
I am mapping various archaeological sites as points; as this is a very rough map, some of the sites share the same coordinates. The points are categorized by the type of site (as "category 1", e.g. settlement, fortifications sanctuary, necropolis etc.). Sometimes, at a single location there are multiple points/datasets per category, like multiple types of fortifications ("category 2").
I chose Point Displacement to declutter the symbols, as seen below. In a next step I just want to show only one symbol per category 1.

In the example below there are three types of fortifications (the tower symbol), like actual towers, walls, fortresses etc. But I just want it to display a single tower-symbol to show that there are fortifications in general (not showing all of them).

Is there any way to do this?


texts: Using "ogr2ogr -sql" without an input fiile

texts:
Lets say I just wanted to turn WKT to GeoJSON on the command line I could do this:

/opt/gdal-custom/bin/ogr2ogr -f CSV /vsistdout/  foo.sqlite \
-dialect sqlite \
-sql "SELECT AsGeoJSON(GeomFromText('POLYGON ((-94.1748 26.90248,-94.1748 27.60567,-93.07617 27.60567,-93.07617 26.90248,-94.1748 26.90248))', 4326)) AS wktIntersection"

However I don't actually need the sqlite file. Is there some sort of null datasource I could use? I tried replacing foo.sqlite with /dev/null and /vsimemory/. Its just a minor annoyance.


texts: How to make dark gray background on TIF image transparent or white in QGIS

texts:
see also: Imported TIF has dark gray background instead of white

That question addresses why the TIF imported into QGIS has a dark gray background when the color set is set to 4096 when exported from PDFCreator.   That file's gdalinfo shows NBITS=4

Here, in this question, I would like to know if there is a way to render the dark gray transparent, much in the way that entering a "0" in
properties->transparency->no data value->additional no data value
will make the black lines transparent.

Alternatively, is there a way to have QGIS render the dark gray as white?


texts: Metrics of change between vector datasets

texts:
I have an old and a new set of parcel data which I'm looking to compare. I'm looking to compare changes in individual parcel geometries, so boundary/corner adjustment. Currently each parcel we have has a unique ID which is updated when it is subdivided or merged, so divisions/mergers are known and there are no other attributes to change.

Are there any standard metrics to measure the changes between each pair of parcels? Ideally I'd like to compile a report over a region and report on a moving average of change.


texts: How to create geometry from set or array of lat/lon in table using PostGIS?

texts:
My question is similar to Creating geometry from lat/lon in table using PostGIS? , but I need to convert set or array of LON & LAT in PostgreSQL into geometry value.

For example, I can convert single LON & LAT to geometry point using below code (using PostgreSQL command) :

SET geometry = ST_GeomFromText('POINT(' || lon || ' ' || lat || ')',4326);

What if I have array of LON & LAT like this lon1,lon2,lon3,lat1,lat2,lat3 and I want to combine those coordinates to geometry values which will produce a long polygon shape instead of a single point shape?


texts: Zoom dependent layers and properties in Leaflet

texts:
I have a GeoJSON file with different values for the "type" attribute and I am trying to open them depending on the zoom. I have tried several possibilities but I have not been able to.

The html is like this:

L.geoJson(river, {
      onEachFeature: function (feature, layer) {
                layer.setText(feature.properties.N, {center: true,offset: -5, attributes: { fill: 'blue' }});
            },

        style: function (feature) {

        var type = feature.properties.T
                if (type == 1) {
                return {
                color: " #ff0000 ",
                weight: 5,
                opacity:1
                }
            }
                if (type == 2) {
                return {
                color: " #088da5 ",
                weight: 5,
                opacity:1
                }
            }

                if (type == 3) {
                return {
                color: " #FFA500 ",
                weight: 3,
                opacity:1
                }
            }

                if (type == 4 ) {
                return {
                color: " #00B200 ",
                weight: 2,
                dashArray: '4, 4',
                opacity:1
                }

                };
           },

    }).addTo(map);


texts: Isolate contour line with a minimum number of points (to eliminate small polygons) with ArcGIS

texts:
I extracted a single contour line (altitude 3660 m) with ArcGIS. How can I eliminate small polygons ?
I know how to do it with a GRASS command (generate contour lines) which simply asks for a minimum of points to create the contour line. I would like to know how do it with ArcGIS.

Here is a red arrow showing the contour line I want to see. The yellow arrow points the small contour lines I want to erase.


texts: Performing Maximum Likelihood Classification in R

texts:
I'm trying to perform maximum likelihood classification of landsat imagery using R, I have been searching for a package that implements it but so far I haven't found it.

What is the easiest way to achieve this? I know package rasclass has some maximum likelihood classification method but what it receives as input is a format completely different than the one I have, as far as I can see.

Is there a package that allows me to do this, or a code snipped anyone could provide? I have been looking at rasclass' source code and thinking about modifying it to do receive the format I want as input, but I'm afraid it's a bit too complicated and I'll make a mistake somewhere.


texts: RiverGIS under QGIS - Not able to execute "Extracting Manning's n values for cross-sections"

texts:
I am using RiverGIS plugin in QGIS and facing a problem in Extracting Manning's n values for cross-sections. The error code is as follow: -

Extracting Manning's n values for cross-sections
UndefinedFunction('function st_line_interpolate_point(geometry, numeric) does not exist\nLINE 39:         ST_Line_Interpolate_Point(inter_xs_dump.geom, 0.0000...\n                 ^\nHINT:  No function matches the given name and argument types. You might need to add explicit type casts.\n')
Process aborted!

I am not very expert in in coding, so please help me.


texts: Change style properties of several layers simultaneously?

texts:
I'm working on bulk data. I have a few hundred GeoTIFF images imported into a QGIS project, and I'd like to modify their layer properties - things like the no-data-value, transparency of certain pixel values, colour map assignments...

But I can't seem to find a way to apply these changes to several layers at once, short of editing the XML file with a global search-and-replace. Is there a way to do this using the GUI?


texts: Extracting raster using values in ArcPy?

texts:
I have an original .tif file which has values ranging from 2 to 28. From that original .tif file, I am extracting cells having a value of >=10 using the following code:

import arcpy, os
from arcpy import env
from arcpy.sa import *

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
env.workspace = "C:/Subhasis/Test/2-25-16/HSA"

outws="C:/Subhasis/Test/2-25-16/HSA/HSA_10"

#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")

inraster = arcpy.ListRasters("*", "TIF")

for i in inraster:
    flds = ("VALUE", "COUNT")
    dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(i, flds)}
    sumcnt = sum(dct.values())
    dct1 = {k:v for (k,v) in dct.items() if k >= 10}
    sumcnt1 = sum(dct1.values())
    percentage=(float(sumcnt1)/float(sumcnt))
    print i,percentage
    newraster = ExtractByAttributes(str(i), "VALUE>=10")
    outname=os.path.join(outws,str(i))
    newraster.save(outname)

The code is working. However, when I open the new extracted tif file in GIS, I found that the cell values range from 10 to 30 instead of 10 to 28.


texts: Simple csv being correupted when imported as Google Earth Engine asset

texts:
I have a CSV of ~500 lat/long locations that I import into my assets as a Table. When I import that Table into my GEE script as a FeatureCollection, the lat/long values get distorted a bit after 5 decimal places. Why is that happening and how do I avoid that?

i.e. the csv contains the lat/long location [7.908736706, 37.2028923] and the imported GEE table contains [7.90873743310952, 37.20289098574329]. Close but no cigar! What's going on here?


texts: How to batch convert Shapefiles to GMLs?

texts: I know how to do this within Quantum GIS but I have 30+ shapefiles that need to be converted to GML.  Is there a Python script sample or plugin for Quantum GIS that I can use to convert my shapefiles to GML?


texts: Importing Google coordinates to QGIS but can't select X/Y field, because coordinates are shown in one row

texts:
I am new to QGIS. I am attempting to import CSV file with x/y coordinates into QGIS. When I add CSV file to delimited text layer all coordinates are shown in one row, so I can't select which one is X or Y. It's the same when I import txt.file.

I am working on Mac so maybe that's a problem. I also saved a file as Windows formatted file CSV. This has not helped ethier. The X and Y coordinates are in columns, but when I have to select X or Y field the columns are merged together as X;Y.


texts: Export a map into a useable format for a design software - QGIS?

texts:
I have a general map of the Mediterranean Sea, showing all the countries I've considered in a study, with watersheds etc. I've tried to make it look nice, but I think it could do with a bit of extra design. All the features in the map are contained in seperate shapefiles (countries, watersheds etc).

Is there a way to export this map, either in a shapefile format or using the print composer export types (PDF, SVG, image) so that it can be enhanced with a design software such as Indesign, Illustrator, or something similar? I'm no designer, just wandering if there is substance for improvement.


texts: Using ArcGIS to split a raster DTM of a river valley into equal parts/sections?

texts:
I would like to split a raster DTM of a river valley into several (maybe 30) equal sections. The raster currently follows the river channel and so is not a straight line but bends and moves like a river. Is there a way for me to split this raster into equal sections which can then be analysed individually? I tried the split and divide tools but these only work for polygons/lines. I am unsure how to divide the raster equally and/or into a specific number of sections. I almost want it to be in a grid-like format.
So for example, I have a section of the river:

And would like to split it up (for example) as so - I did this by calculating the Euclidean Distance:

As if each colour was its own section that I can analyse individually.


texts: Running SAGA Processing algorithm from console takes much much longer than from GUI?

texts:
I have a road layer with 32000 features, and I want to find all the intersections. when I use processing->SAGA -> convert polygon/line vertices to points, it is all done in less than  5 seconds. However, when I try to do the same thing in PyQgis:

vlayer = QgsVectorLayer("Point", "temp_points", "memory")
c=processing.runalg("saga:convertpolygonlineverticestopoints",layer,vlayer) # get intersection points

it takes about 15 minutes.


texts: Unable to export data to ShapeFile in ArcGIS

texts:
I've imported an Excel sheet into my ArcGIS project but am unable to convert this to a ShapeFile. When I try, via export data, I get an error message saying either 'output name is not valid', whatever I set the name as, or that it simply cannot save the file.

Any ideas on where I am going wrong?


texts: QGIS creating user CRS rather than using UK National grid

texts:
Does anyone have any issues using QGIS (current version 2.12.1 on mac) with the British National Grid CRS (EPSG 27700)?

All of my data is held in BNG and my workspaces are set up to default to BNG when I create new datasets. However, almost every time I create a new dataset (raster or vector) QGIS defaults to creating a user generated CRS; an example in the image below.  These are always slightly out of sync with the 27700 data but shift back if you change the CRS to 27700.  Sometimes I'll finish working with 10 different user CRSs.

Is this something that I'm doing or is this a bug?  I can't remember whether this has ever happened on previous versions of QGIS.


texts: ArcMap correct Errors in Attribute Table

texts:
I want to replace the all the -9999 values in my attribute table with 0.

Does anyone know an easy way to do this?


texts: Irrigated Area Recognition in Google Earth Engine

texts:
I'm using Google Earth Engine with the Python API.
I have calculated the NDVI value over an area, set a threshold and the result is a layer where black areas exceed my NDVI threshold (don't pay attention to the red and blue dots, they are not important in this issue)

Now, I want to automatically identify these black areas (the ones that exceed a given extension) and save their position in order to get a list of the (probably) irrigated camps in this ROI.

How could I do this?


texts:
I need to filter the WMS layer included in my map for a date field. Basically I only want to show what has happened after a certain date:

I tried the following, but no luck

createMap = () => {

    addCommonProjections();
    const vectorSource = new VectorSource();

    const vectorLayer = new VectorLayer({
        source: vectorSource
    });

    const baseLayer = new TileLayer({
        source: new OSM()
    });

    this.wms_source = new ImageWMS({
        url: environment.live_lightning_data,
        projection: 'EPSG:4326',
        styles: 'default',
        params: {
            layers: '0'
        },
    });

    const lightning_layer = new ImageLayer({});
    this.map = new Map({
        layers: [baseLayer, lightning_layer],
        target: document.getElementById('lightning'),
        view: new View({
            center: fromLonLat([175.79, -37.79]),
            zoom: 7,
            projection: 'EPSG:3857'
        })
    });

    lightning_layer.filter.greaterThan('data_date_time', '2019/03/02');
    lightning_layer.setSource(this.wms_source);

I know, a little chaotic, but the map map itself displays the data correctly, but it is unhappy with however I try to add the filer.

The way I added it now,results in '... is not a function'

Anyone any suggestions?


texts:
Is there a way in FME to process only, say, 100 records in a dataset? While developing a workspace I often have to iterate until I get it right. If the source dataset is large it's a waste of time to watch it process the whole thing just to do it again after tweaking a setting.

I'm thinking of something like using the counter transformer and a test, but that would still (partially) process all of the input data.


texts: Insert several gtfs feeds into the same postgis/pgrouting database

texts:
I want to import sevaral gtfs feeds into my postgis/pgrouting database so that I can route over all my gtfs feeds.

The problem is that the sources / destinations ids of my gtfs are likely to conflict.

How can I solve this conflict issue?


texts: Rotating and moving polygons to different angles using ModelBuilder/ArcPy?

texts:
I have a few thousand polygons (all square) that I need to rotate to various angles that are in the attribute table, and then move so that each polygon gets moved 3/4 up from the original centroid, which I have as a point shapefile. It would look like this:

Ideally, I would like to build a model for this process so I can repeat it in the future. I am not all that comfortable with Python, but if that is the only way to get the desired result, I will be willing to try it. I have arcGIS 10 with an ArcView license.


texts: Draw icons and lables only for certain entities in Tilemill 2

texts:
I've recently started to use mapbox' TileMill2 to play around with vector tiles. It's pretty nice and, in comparison to Tilemill(1), much more faster. Exactly what I was searching for.

Now I'm trying to visualize certain entities with their icon and label, for example only police stations and restaurants. Placing the icons only for those features works, but not the labeling part. Although the labels for the chosen entitites are placed, but unfortunately also all other. Maybe you have an idea on how to solve my problem? I've already checked the references, but no success.

Find below the corresponding css stuff and a picture )right click on picture and "open in new tab" for bigger pic), which shows the problem:
[-] Only restaurants and police stations got visualized by icons (this is exactly the behavious, which I want to have)
[-] Restaurants and police stations got labeled. This is also the right behaviour.
[-] Everything else is labeled. I don't want this :-/ I've thougt that the ::label-part would solve this. But probably there is an error.!

#poi_label[zoom>=15] {
  // Separate icon and label attachments are created to ensure that
  // all icon placement happens first, then labels are placed only
  // if there is still room.

  ::icon[maki='police'] {
    // The [maki] field values match a subset of Maki icon names, so we
    // can use that in our url expression.
    // Not all POIs have a Maki icon assigned, so we limit this section
    // to those that do. See also <https://www.mapbox.com/maki/>
    marker-fill:#666;
    marker-file:url('icon/police-12.svg');
  }

  ::label {
    text-name: @name;
    text-face-name: @sans_md;
    text-size: 10;
    text-fill: #f55;
    // text-transform: uppercase;
    // text-character-spacing:  0.25;
    // POI labels with an icon need to be offset:
    [maki='police'] { text-dy: 8; }
  }

  ::icon[maki='restaurant'] {
    // The [maki] field values match a subset of Maki icon names, so we
    // can use that in our url expression.
    // Not all POIs have a Maki icon assigned, so we limit this section
    // to those that do. See also <https://www.mapbox.com/maki/>
    marker-fill:#666;
    marker-file:url('icon/restaurant-12.svg');
  }

    ::label {
    text-name: @name;
    text-face-name: @sans_md;
    text-size: 10;
    text-fill: #f55;
    // text-transform: uppercase;
    // text-character-spacing:  0.25;
    // POI labels with an icon need to be offset:
    [maki='restaurant'] { text-dy: 8; }
  }
}


texts: Determining degrees of freedom for Geographically Weighted Regression (GWR) in ArcMap?

texts:
I'm doing Geographically Weighted Regression (GWR) analysis in ArcMap and need the degrees of freedom for the model in order to calculate p values.

The summary does not list degrees of freedom, but lists an "effective number".

Is this the same thing?


texts: What are pros and cons of Python add-ins vs .NET add-ins to ArcGIS for Desktop?

texts:
I found in ESRI .NET help:

Python
All ArcGIS Desktop applications include an embedded scripting language called Python. Many areas of ArcGIS—particularly geoprocessing—are accessible through simplified Python application programming interfaces (APIs), making it easy to author and automate common tasks. Python scripts are easily shared and can be produced without an external development environment. There are a variety of public domain Python modules focused on areas such as science, engineering, and mathematics.
With all its strengths, Python is not suited for every programming task. Although coverage is improving, not all areas of ArcGIS are currently exposed to Python. In addition, the editing and debugging experience in Python is not as sophisticated or easy to use as those included with commercial development environments such as Visual Studio. Lastly, you cannot listen for and respond to ArcGIS events, implement a COM interface, or plug into ESRI’s many COM extensibility points using Python.

and in ArcGIS Desktop help:

ArcGIS 10.1 introduces Python to the list of languages for authoring Desktop add-ins, providing you with an easy solution to extend desktop functionality. To simplify the development of Python add-ins, you must download and use the Python Add-In Wizard to declare the type of customization. The wizard will generate all the required files necessary for the add-in to work.

If Python and .NET add-ins have the same functionality?
Does that mean that I can do the same thing in Python  as .NET ??


texts:
I have found out how to get tab-delimited X,Y data into a layer but the Z field value is not given as an option on import?

I tried to edit the attributes to add and define a new field but it keeps crashing.


texts:
According to web posts from a few years ago, exporting to .eps was easy with QGIS.  Just choose : Export as image --> select .eps type.  However, I do not find that functionality in the current version of QGIS.
Has it been removed as an option or is there another route to get there?
Context: I am preparing my figures for submission to a journal for publication.  The journal requires over 300dpi (high resolution) .tiff or .eps format.  When I export my map from QGIS in over 300dpi as a .tiff, the file size is too big for the journal's webserver to upload.  As an alternative, I am hoping to export as .eps to conserve file resolution but keep file size down as I have with my graphs exports from R.


texts: Error: "ee.ImageCollection(...).subtract is not a function" (Earth Engine)

texts:
When I execute this expression:

var ETfraction = (Thot-LST)/(Thot-Tcold)  // **pseudo code**

I get an error:

  ee.ImageCollection(...).subtract is not a function

Why does this not work? Here is my code:

//Find Reference ET (ETo) from Terraclimate

var dataset = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
              .filter(ee.Filter.date('2015-01-01', '2017-12-30'))

.map(function(image)

{return image.clip(basin)});

var pet  = dataset.select('pet');

var pet = pet.map(function(img) {

 return img
    .multiply(0.1)
    //.subtract(273.15)

.copyProperties(img, ['system:time_start']);
});

var vis = {

min: 50,

max: 200.0,

 palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

Map.addLayer(pet,vis);

print(pet,'PET')

//Find LST from MODIS(MOD11A1.006 = Terar land surface temperature daily 1 km.)

var dataset = ee.ImageCollection('MODIS/006/MOD11A1')

.filter(ee.Filter.date('2015-01-01', '2017-12-30'))

.filterBounds(basin)

 .map(function(image){return image.clip(basin)});

var lst = dataset.select('LST_Day_1km');

var lst = lst.map(function(img) {

 return img

 .multiply(0.02)
    // .subtract(273.15)

.copyProperties(img, ['system:time_start']);
});

var landSurfaceTemperatureVis = {

 min: 14000.0,
  max: 16000.0,

 palette: [

  'blue', 'limegreen', 'yellow', 'darkorange', 'red'
  ],
};

//Ts=LST of each pixel

var monthly_lst = ee.List.sequence(0, 3*12)

.map(function(n)
{

// .sequence: number of years from starting year to present

var start = ee.Date('2015-01-01').advance(n, 'month'); // Starting date

var end = start.advance(1, 'month'); // Step by each iteration

return lst.filterDate(start, end)

.mean()

 .set('system:time_start', start.millis());

});

var collection = ee.ImageCollection(monthly_lst);

//Thot = LST for the hot pixel

var Th_lst = ee.List.sequence(0, 3*12)
.map(function(n) {

// .sequence: number of years from starting year to present

var start = ee.Date('2015-01-01').advance(n, 'month'); // Starting date

var end = start.advance(1, 'month'); // Step by each iteration

return lst.filterDate(start, end)

        .max()

        .set('system:time_start', start.millis());

});
print(Th_lst);

var collection2 = ee.ImageCollection(Th_lst);

 print(collection2,'Thot');

 Map.addLayer(collection2,vis);

//Tcold = LST for the cold pixel
var Tc_lst = ee.List.sequence(0, 3*12)
.map(function(n) {

// .sequence: number of years from starting year to present

var start = ee.Date('2015-01-01').advance(n, 'month'); // Starting date

var end = start.advance(1, 'month'); // Step by each iteration

return lst.filterDate(start, end)

        .min()

        .set('system:time_start', start.millis());

});

print(Tc_lst);

var collection3 = ee.ImageCollection(Tc_lst);

 print(collection3,'Tcold');

 Map.addLayer(collection3,vis);

// // //ET fractional (ETf=(Thot-LST)/(Thot-Tcold))

var etf = ee.ImageCollection(Th_lst)

          .subtract(ee.Image(monthly_lst))

          .divide((ee.Image(Th_lst))

          .subtract(ee.Image(Tc_lst)))

         .rename('ETF').cast({'ETF': 'float'});

print(etf);

Map.addLayer(etf,vis);

Code Editor link


texts:
I have a polygon layer and I want to know the altitude of the centroids. I guess I have to find on the internet (how?) a layer that contains data about the territory's height (in this case Italy in particular Toscana and Liguria) and to find a way to correlate it with my layer (using Python is better). Does anyone have a suggestion?

vlayer = QgsProject.instance().mapLayersByName('Centroids')[0]
rlayer = QgsProject.instance().mapLayersByName('Elevation')[0]
params = {'INPUT_VECTOR': QgsProcessingFeatureSourceDefinition(vlayer.id(), True),
          'INPUT_RASTER': QgsProcessingFeatureSourceDefinition(rlayer.id(), True),
          'COLUMN_PREFIX' : 'Altitude',
          'OUTPUT': 'memory:'}

newLayer = processing.run("qgis:rastersampling", params)
layerOutput = newLayer['OUTPUT']
QgsProject.instance().addMapLayer(layerOutput)

I tried in this way, but there's something wrong. This is the error message: Could not load source layer for INPUT: no value specified for parameter


texts: Determining if layer is visible using PyQGIS

texts:
I am working on a plugin 'Stripchartdrawer' current commit https://github.com/sickel/qgisstripchart/tree/ba3afd7184e33331b5af38111e91dbb9b6a2d1c0.

One of the features in this plugin is to be able to visually select features in a layer based on the value of a field by clicking or dragging in the chart. There are two cases when this kind of selection does not work. One is if the layer does not have a unique id field, if so, I am able to give a warning.

The other case is if the layer is not activated, i.e. checked in the layer panel, I would like to give a warning in that case as well if it seems as the user tries to select data, but how can I find out if the layer is activated?

In my mousereleaseevent I am doing

self.layer.select(selectedids)

Where selectedids is one or more ids that I am to select.

What I want to achieve is something like (pseudocode:)

if self.layer.isActive():
   self.layer.select(selectedids)
else:
   self.iface.messagebar.pushmessage('Warning',(...))

But what is the real "isActive()"?

I am not 100% sure about the nomenclature here. a layer can be not visible either because it is not checked in the layer panel, if the view is outside the extent of the layer or if the layer is covered with another layer, one layer at a time can be activated, ie eg being the layer "Select features by plygon" will operate on. I want to be able to check if a layer is checked (no pun intended) and do not care if it is activated as described above.


texts: What can I do with lidar returns of 6 or greater in LAS 1.2 specification?

texts:
I have lidar LAS data that contains a number of points that represent returns 6 or 7.  The data is supposed to be in ASPRS LAS 1.2 specification, which does not allow return numbers (echos) greater than 5.  With returns of 6 and 7, my data is not 1.2 spec compliant.

For all of you point cloud junkies, what are my options to get the data into 1.2 compliance, without deteriorating its integrity?

EDIT:

LiDAR Sensor: Riegl LMS_Q680i

Recording Software: Riegl Ri-Acquire

Raw Data Processing Software: Riegl RiANALYZE 6.0.2 (build: 2013-12-04)

Screen cap showing header of sample data with lasinfo


texts: pyqgis correct activate signal for combo box list

texts:
I am a pyqgis user and I am trying to create a QGIS plugin where that plugin have two tabs widgets and two combo box where first combo box have a list from layers and second combo box have list of fields for current layer in combo box one.

That combo box is inside in tab widget two, and I want to to activate any time user start plugin.

My code works but not correctly - I have curious results.

In the first way in my _ init _  the combo box 1 work fine every time to start the user the plugin then combo box 1 activate correct, but the combo box 2 insert only six fields for some layers (that layers have many fields 19).

In the second way combo box 1 and combo box 2 work fine have correct list of layers and fields, but to activate combo box 1 must be the user to move some layers in panel and after this activate. I don't want this .

Where I have wrong ?

Who is correct signal activate for combo box 1 ?

Maybe some signal to activate when the user press tab 2 from tab widget?

in _ _init __  :

#first way
registry = QgsMapLayerRegistry.instance()
registry.layersAdded.connect(self.change_layers)
#second way
self.iface.legendInterface().groupRelationsChanged.connect(self.change_layers)

#field change
self.ui.comboBox_2.currentIndexChanged.connect(self.updateFields)

layer function :

def change_layers(self):
   self.ui.comboBox_2.clear()
   legendInterface = self.iface.legendInterface()
   listLayerName = [i.name() for i in legendInterface.layers() if i.type() == QgsMapLayer.VectorLayer]
   # add all these layer names to combo box comboLayerList
   self.ui.comboBox_2.addItems(listLayerName)

field function :

    def layer_field(self):
        self.mapCanvas.refresh()
        layers = self.iface.legendInterface().layers()
        # Identify selected layer by its index
        selectedLayerIndex = self.ui.comboBox_2.currentIndex()
        print selectedLayerIndex
        selectedLayer = layers[selectedLayerIndex]
        # Identify fields of the selected layer
        fields = selectedLayer.pendingFields()
        # Get field names of the fields
        fieldnames = [field.name() for field in fields]
        # Clear comboBox_5
        self.ui.comboBox_3.clear()
        # Add field names to comboBox_5
        self.ui.comboBox_3.addItems(fieldnames)


texts: My self-hosted CartoDB map references /api/v1/map -- No route matches

texts:
I'm using a self-hosted install of CartoDB. When I try to publish a visualization using the sample Europe data the visualization page references /api/v1/map which results in a routing error.

ActionController::RoutingError (No route matches [GET] "/api/v1/map"):

I'm serving CartoDB with Apache/Passenger, and all of the other endpoints seem to be working.

I have a ProxyPass setting in Apache to proxy /api/v1/sql requests to the CartoDB SQL API service, but /api/v1/map (and other /api/v1 requests) are passed through to CartoDB itself.

I'm not sure what other info would be useful. Please let me know and I'll provide it.


texts:
I have a point layer with multiple points and a layer with the starting point (fire station).

In the end, I want to create a contour file with distances. But in such a way that the approximate distance is given and not the direct distance.

My question is:
How can I create the lines between points and the lines must be connected to the nearest point.

I want wo be able to measure the distance approximatley by road and not the direct distance (see below).


texts: Is there any way to get an email when ArcGIS Server services are stopped/down?

texts:
Occasionally, but all too often unfortunately, my AGS services just stop for no reason. I don't usually notice this until the users call me and ask why their maps have no data.

Is there any way (I was thinking Python) that I can get a response from AGS when services stop? How do other AGS users deal with this problem?


texts: How to extract sublinestring from a linestring with overlapping segments

texts:
Problem statement:
Given a bus route (defined as a linestring ln.the_geom), and a set of bus stop locations (pt1.the_geom,pt2.the_geom,...ptn.the_geom), divide the bus route into segments between bus stops.
I am using two postgis functions

    ST_line_substring(ln.the_geom,ST_line_locate_point
    ln.the_geom,pt1.the_geom),ST_line_locate_point(ln.the_geom,pt2.the_geom))

The above works well if the bus route has no overlapping sections. However, if there is overlapping sections, the ST_line_locate_point will snap a point to the other end of the linestring.
As an example: Bus travels from A to B, make a loop, then travels from B to A.
"A to B" and "B to A" are parallel routes with all points overlapped. In case of bus stops P1 and P2, P2 is mapped into point P2' instead (on the opposite direction). The segment is extracted wrongly in this case (where extra segments are included).

How to solve problems like this? The overlapping bus routes can happen at any section of the route...


texts: Error 010240 while running cell statistics operations

texts:
Can somebody figure out whats going on here?
I had absolutely no problems running these commands before.
I am trying to add 1046 raster layers. I created a list of rasters using arcpy.ListRasters()
Then used CellStatistics to "SUM" all the rasters in the list. Error 010240 was generated. Then repeated the procedure for only a few rasters. It worked. Repeated again for all the rasters and the error was generated again. Came back to a few rasters and the error is generated again! Clicking on the error either closes Arc with 'Not enough memory' message or nothing happens. Searching a tool and opening it opens a blank 'Script error' message. I'm using Arc 10.3. I recently installed Point cloud library and CMake and all the other dependencies. Both simple addition by looping and cell statistics are not working in any case now.

>>> ws="E:\Karun\MTech_Project\MB_RasterAnalysis.gdb"
>>> arcpy.env.workspace=arcpy.env.scratchWorkspace=ws
>>> rlist1
    [u'rast_1', u'rast_2', u'rast_3', u'rast_4', u'rast_5']
>>> fr=CellStatistics(rlist1,"SUM")
    RuntimeError: ERROR 010240: Could not save raster dataset to E:\Karun\MTech_Project\MB_RasterAnalysis.gdb\sumco_ras with output format FGDBR.
>>> fr=CellStatistics(rlist1,"SUM")
    RuntimeError: ERROR 010240: Could not save raster dataset to E:\Karun\MTech_Project\MB_RasterAnalysis.gdb\sumco_ras with output format FGDBR.
>>> rlist1[0]+rlist1[1]
    u'rast_1rast_2'
>>> Raster(rlist1[0])+Raster(rlist1[1])
>>> temp=Raster(rlist1[0])+Raster(rlist1[1])
    RuntimeError: ERROR 010240: Could not save raster dataset to E:\Karun\MTech_Project\MB_RasterAnalysis.gdb\plus_ras with output format FGDBR.
>>>


texts: Adding new record with specific coordinates into shapefile using ArcGIS for Desktop?

texts:
I need to create a new point record in my shapefile. When I use the Create Features tool in the Editor tool, I can click anywhere on my map to create a point, but when I go into the attribute table to change the point to exact lat/long I want, the point does not move to those coordinates.

Any ideas on how to add a point record into my shapefile with a specific lat/long properly?


texts: Whitebox's cost distance tool stuck in endlesss loop

texts:
Cost distance function from WhiteboxTools (v.0.11.0) gets stuck in an endless loop after execution. The inputs are 8-bit unsigned interger GeoTIFFs. The source raster has values 1 and Null (size: 16 KB), and cost raster has values from 10-54  (size: 460 KB) (see below for an image of source (grey), cost raster (gradient white-to-red) and white is null/nodata). If necessary, I can provide the GeoTIFFs if needed.

WhiteboxTools - Cost distance code:

import sys
wbtools_dir = r'C:\WhiteboxTools_win_amd64\WBT'
sys.path.insert(0, wbtools_dir)

from whitebox_tools import WhiteboxTools

wbt = WhiteboxTools()

source = "c:/src2.tif"
cost = "c:/cst1.tif"
out_accum = "c:/cd.tif"
out_backlink = "c:/backlink.tif"

wbt.cost_distance(source, cost, out_accum, out_backlink)

Output:

***************************
* Welcome to CostDistance *
***************************
Reading source data...
Reading cost data...
Initializing: 0%
Initializing: 1%
...
Initializing: 99%
Initializing: 100%
Loop 1: 0%
Loop 1: 1%
Loop 1: 2%
Loop 1: 3%
...
Loop 4000: 1%

It takes 2-3 seconds to complete each loop and the process was killed after 30-45 minutes. I replicated the cost distance in ArcMap and run time was ~45 seconds.

The code comes directly from the WhiteboxTools user manual and it executes without any issues, so I suspect it is either user error (need to format my source and cost rasters in certain way) or Whitebox bug. Any ideas? Also, I am able to execute other Whitebox tools such as breach_depressions() and create valid outputs.


texts: Gdal Python Bindings and File Geodatabase API

texts:
I'm trying to get python 2.7 to work with GDAL/OGR and an ESRI file geodatabase.  There have been a few similar questions asked but the answers don't fully resolve my situation.

I have installed the 64-bit GDAL Core files from http://www.gisinternals.com/sdk/ I have also installed the Python Bindings and the File GDB.

The background to the problem is that when I run ogrinfo --formats from the command line it returns: "FileGDB" (read/write), when I run ogr from the command line I can see the layers in the Geodatabase BUT when I try to use the following code in Python in returns "None"

driver = gdal.GetDriverByName('FileGDB')
ds = driver.Open(r"my file.gdb",0)
print driver

Ogr definitely works in python because the following code returns the correct response:

driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(r"My Shapefile.shp",0)
print ds

driver = gdal.GetDriverByName('GTiff')
ds = driver.Open(r"My GeoTiff.shp",0)
print driver

Does anybody have any pointers to resolve this.

Thanks


texts: What unit is the Shape_Length field in for the GCS_North_American_1983_CSRS projection?

texts: How can I convert these to values to meters? I'm trying to use the Densify tool but can't because I'm assuming the distance field is tied to Shape_Length.


texts: Adding nearest point-(features) to attribute table

texts:
Using QGIS I'd like to add the (10) nearest random points, that I created within a polygon, to a target point and have their distances to said point added to that point's attribute table.
I'd like to do this for multiple target points, but cannot figure out how to do it.

this is how I'd like my attribute table to look like (used Excel only to visualise).


texts: QGIS 2.8 & Python: Merge problems using Python

texts:
I have 9 vector layers that I rasterized. All of these have the same:

x;y origin
pixel size (1 pixel=1 meter)
SRS

When displayed separately in QGIS these 9 rasters are superimposed correctly but when merge using the Python script there is one or more rasters that are offset approximately 100 meters to the North.

Can you check if the script is correct or if the problem come from my data?

import os
import gdal
import struct
import numpy
import osr
import sys

def raster2array(raster):
    band = raster.GetRasterBand(1)
    return band.ReadAsArray()

def getOrigin(raster):
    geotransform = raster.GetGeoTransform()
#    print geotransform
    return [geotransform[0], geotransform[3]]

def getLayers():
    return iface.legendInterface().layers()

def array2raster(rasterfn,newRasterfn,array, xPadding = 0, yPadding = 0):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols + xPadding, rows + yPadding, 1, gdal.GDT_Float64)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

gdal.UseExceptions()

layers = getLayers()
firstLayer = layers.pop()
path = firstLayer.dataProvider().dataSourceUri()
raster = gdal.Open(path)

print "Using \"" + firstLayer.name() + "\" as base layer"
origin = getOrigin(raster)

dest = raster2array(raster)

for layer in reversed(layers):
    print "Processing \"" + layer.name() + "\"."
    provider = layer.dataProvider()
    path = provider.dataSourceUri()
    raster = gdal.Open(path)
    layerOrigin = getOrigin(raster)
    layerData = raster2array(raster)

    for x in range(0, len(layerData)):
        for y in range(0, len(layerData[x])):
            if layerData[x][y] > 0:
                dest[x][y] = layerData[x][y]

rasterfn = '/Users/sbooob/Desktop/QGis/SIG/1CLC12_Raster.img'
newValue = 0
newRasterfn = '/Users/sbooob/Desktop/QGis/SIG/out.tif'

array2raster(rasterfn, newRasterfn, dest)


texts: QGIS 2.2 - possible to do a JOINS from the Action list?

texts:
I have a set of layers that need to be joined (Layer Properties > Joins) to one specific layer. The reason for this is that I will be updating and adding new layers into my project and would like a single Action command for the specific layer to automatically have the selected layers join with it. I want to avoid having to run joining/merging algorithms altogether.

Is this easily achievable?

After spending time searching on how to create Actions, most of what I found describe the same command which is how to open up a web browser to search for whatever is in the specified attribute. Unfortunately, the default actions don't seem to contain anything similar to what I am looking for.