texts: Processing raster files with raster calculator in PyQGIS?

texts:
I am trying to process some raster files using python scripts in Qgis 2.20, and I have been having trouble finding any documentation on using raster calculator in scripts being run in the processing toolbox.

I am trying to use scripts to process weather data. I will need to recalculate values of cells in order to represent projected wind speeds at different heights and probable gusting speed. I am trying to use the raster calculator for this but have been having trouble. I was not able to easily find a function for the script, and after some digging through the python console I was able to find the saga:rastercalculator function.
Documentation in the python console (obtained through alhelp()), say that the parameters are as follows.
algorithm: raster calculator
Grids
formula
Raster

I am trying to run a basic formula through the raster calculator to get it working, but am having no luck.

Is this even the right approach to recalculating grid values in a python script?

Here is the code:

#load and symbolize tiff

##[User scripts]=group
##input_raster=raster

#Imports
#=======
from qgis.core import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.analysis import  *
import processing
from processing.core.Processing import Processing

#import Tiff

fileName = input_raster
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
raster = QgsRasterLayer(fileName, baseName)
raster.loadNamedStyle(r"C:\Users\Maximilian\Documents\School\AGRG\test.qml")
# add map to display
QgsMapLayerRegistry.instance().addMapLayer(raster)

#recalculate values
processing.Saga:RasterCalculator(raster, "(raster)*2/2",outraster)

Here is the error:

Uncaught error while executing algorithm
Traceback (most recent call last):
Traceback (most recent call last):
  File "C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins\processing  \core\GeoAlgorithm.py", line 202, in execute
    self.processAlgorithm(progress)
  File "C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins\processing\script\ScriptAlgorithm.py", line 252, in processAlgorithm
    exec script in ns
  File "<string>", line 65
    processing.Saga:RasterCalculator(raster, "(raster)*2/2",outraster)
               ^
SyntaxError: invalid syntax


texts: QGIS Points in Polygon Crashes (with data)

texts:
New to QGIS, using 2.14.0 on Windows 7

I am trying to run Vector -> Analysis Tools -> Points in Polygon... on these data. I want to get the PNTCNT and summary statistics on the Amount attribute. However, every time I run Points in Polygon..., I get a "crash dumped" error after reaching 99%:

The only other relevant question I've seen regarding this issue didn't have much in the way of responses. It did have a question in response asking if each layer was a shapefile. In my instance I am trying Points in Polygon... with a shapefile and a .csv layer. I'm not sure if that may be my problem. However, I tried saving my delimited layer as a shapefile and trying this again, with no luck.

Also, I've tried Vector -> Data Management Tools -> Join Attributes by Location... and I get a Python error (screenshot with data).

This all leads me to wonder if there is an issue in my data.


texts: Python 2.7 cannot open CSV I just created....troubleshooting from other posts did not work

texts:
I'm using python 2.7 on windows platform. (I don't have a choice to upgrade - cannot import any non standard modules...work constraints.) I have written a csv file from a feature class. I then am attempting to re-format text in a column. For some reason I am unable to open the csv.
I get the following errors:
D:\JJ_Development\Projects\Basemap\Transportation>new.py
D:\JJ_Development\Projects\Basemap\Transportation
Traceback (most recent call last):
  File "D:\JJ_Development\Projects\Basemap\Transportation\new.py", line 199, in <module>
    newlabel()
  File "D:\JJ_Development\Projects\Basemap\Transportation\new.py", line 97, in newlabel
    with open(rc_existing, 'rb') as csvfile:
IOError: [Errno 2] No such file or directory: 'D:\\JJ_Development\\Projects\\Basemap\\Transportation\\roadcore'

Both the csv and the python script are in my current working directory.
The csv is not read only.
I can open the csv in both notepad++ and arcmap with no issues.
The file is not locked...I have restarted computer and have changed the name as well as copying over a fresh script.

I'm posting a snippet of csv creations (works) and then errors on the line "with open" (obviously this is not the full code)
    #********DATA AND PATHS********
curr_dir = os.getcwd()
print(curr_dir)
arcpy.env.workspace = (curr_dir)

#mxd = arcpy.mapping.MapDocument(r'D:\JJ_Development\Projects\Basemap\Transportation\Labels.mxd')
#in_fc = arcpy.mapping.ListLayers(mxd, 'RoadCore_test2')[0]
fc = (r'D:\JJ_Development\Projects\Basemap\Transportation\Roads.gdb\RoadCore_Existing')

fields = ["ID", "ADMIN_ORG", "RTE_CN", "GIS_MILES"]
identical = ["ID", "ADMIN_ORG", "RTE_CN"]
labels = []

rd_names = (' ROAD', ' STREET', ' DRIVE', ' BLUFF', ' COVE',' LOOP', ' LANE', ' RD', ' LN', ' DR')

#********Initial FS Roadcore Existing Data Cleanup********
rc_existing = os.path.join(curr_dir, 'roadcore')

def test():
    fields = ["ID", "ADMIN_ORG", "RTE_CN", "GIS_MILES"]

    with open('{}.csv'.format(rc_existing), 'wb') as outf:
        dw = csv.DictWriter(outf, fieldnames=fields)
        dw.writeheader()
        with arcpy.da.SearchCursor(fc, fields) as rows:
            for row in rows:
                if row[3] == 0:
                    pass
                elif row[3] == None:
                    pass
                else:
                    dw.writerow(dict(zip(fields, row)))
    print "File written: {}".format(rc_existing)
    outf.close()

#********FORMATTING LABELS********
def newlabel():
    with open(rc_existing, 'rb') as csvfile:
        reader = csv.reader(csvfile, delimiter=',', quotechar='|')
        for row in reader:
            print('csv open')

    # with arcpy.da.SearchCursor(rc_existing, fields) as sc:
        # for row in sc:
            # if (row[0].startswith(' ')):
                # labels.append('{}'.format(row[0].replace(' ', '')))

if __name__ == '__main__':
    test()
    newlabel()


texts:
I was learning eCoginition by tutorial, but still, not familiar with some advanced funcionality, such as stretching raster data.

As you can see, the values are condensed in a small portion of the range. So,the goal is to stretch it to the scale of 0-255 using eCog. And the stretching process should make actual change to the data, not just for visualization.
Any ideas? Thanks in  advance.


texts: Cost Path as Polyline producing empty output using ArcGIS Pro

texts:
I have been trying to execute the Cost Path as Polyline tool using both ArcMap 10.8 and ArcGIS Pro 2.5 with spatial analyst extension.  I'm following the steps in this Esri blog to use Euclidean Distance with barriers to create the distance raster and the out-back direction raster.
The Euclidean Distance tool outputs the distance and out-back rasters as expected.  However, when I run the Cost Path as Polyline tool, the tool executes successfully but with zero geometry.  According to this post, I'm not the first to run into this issue.  Why does the Cost Path as Polyline tool produce empty output?
Distance Raster

Out-Back Raster


texts: Importing Excel of multiple sheets into QGIS

texts:
I have an excel document of multiple separate sheets.

I can import them one at a time by saving the individual sheets as CSV’s but is there away to bring the entire document into QGIS?

The information is point geometries.


texts: Capture non-attribute user input in ArcGIS Pro task

texts:
I am trying to write an ArcGIS Pro task that involves one or more data layers.  However, the exact number of layers may vary between runs. Sometimes the data may be in one file, other times there may be two or more files.

Is it possible to get input from the user that isn't an attribute value?

Specifically, I would like the user to be able to specify with an integer how many different layer(s) are involved, and then use that input value in a For loop containing Add Data command(s).


texts: Converting Point Feature to Latitude/Longitude Coordinate and GPX?

texts: I drew a point feature using the draw tool in ArcMap on a map in an area of interest. I am now interested in extracting the actual latitude/longitude of this point in geographic space so I can actually navigate to this exact area in a car. In addition, I am interested in determining how to export this as a GPX feature and inputting it directly into my Garmin GPS device.


texts: Python- use geocode function on table of addresses & copy lat/long to appropriate row

texts:
I have a table of addresses, and I need my Python script to use my Google API geolocate function to return lat/long coordinates for each and add to a new field in the same row for each address. The geocode function works fine- I just can't get the script to iterate through each row of the table, add the address to the function, and then copy the output lat/long to the field in the same row. here's what I have:

import urllib, json, time, arcpy

arcpy.env.workspace = "D:/GIS/addr.dbf"

#sets variables, adds new field to hold lat/long coordinates
fc = 'addr.dbf'
field1 = 'address'
field2 = 'loc'
arcpy.AddField_management(fc, field2, "TEXT")

#function that uses Google API to geolocate- this part works consistently
def geolocate(address,
api="key_here",delay=4):
  base = r"https://maps.googleapis.com/maps/api/geocode/json?"
  addP = "address=" + address.replace(" ","+")
  gUrl = base + addP + "&key=" + api
  response = urllib.urlopen(gUrl)
  jres = response.read()
  jData = json.loads(jres)
  if jData['status'] == 'OK':
    resu = jData['results'][0]
    finList = [resu['formatted_address'],resu['geometry']['location']
    ['lat'],resu['geometry']['location']['lng']]
  else:
    finList = [None,None,None]
  time.sleep(delay)
  return finList

#adds address field as text to geolocate in function, adds output lat/long
#(indexed locations 0 and 1 from the finList output)
##this is the part that doesn't work!
geo = geolocate(address = field1)
cursor = arcpy.UpdateCursor(fc, [field1, field2])
for row in cursor:
  field2 = geo[0], geo[1]
  cursor.updateRow(row);


texts: How to load features into empty feature class using stand-alone arcpy script

texts:
I have an empty feature class that I would like to populate with features from source feature classes.  I have done this manually using the Object Loader.  However, I would like to automate the process with my stand-alone python script.

I considered using append, but my target and source feature class schemas are different.  I have seen some script for Load Data, however, it is unclear if and how I can create a cross reference database with arcpy.

So, I am wondering if there is an arcpy equivalent to the manual Object Loader.


texts:
Can I create a legend of a map that I had get from a WMS?
Can I create a legend of a raster map with Qgis version 1.7.4, or only with the 1.8 version?


texts: How do I append a shapefile to PostGIS table with multiple geometry fields using ogr2ogr?

texts:
I am having an issue when trying to append shapefile data to an existing table in my PostGIS database using ogr2ogr. The existing table has multiple geometry fields. When I append the shapefile the attribute data uploads fine but the geometry fields are all NULL. I have tested out the upload appending to a new table with only one geometry field and it works fine, but when there are multiple fields it just fails to input any geometry information.
My command looks something like:
ogr2ogr.exe -update -append \
  -f PostgreSQL "PG:host=xxxx port=5432 dbname=yyyy password=***** user=zzz" "C:\xxx.shp" \
  -nlt PROMOTE_TO_MULTI \
  -lco GEOMETRY_NAME=geom \
  -nln schema.tablename \
  -a_srs EPSG:3857

How do I make sure my geometry information is also appended?
I have tried using the sql argument to select columns and rename the geometry column to match the one I want the data in, but nothing has worked so far. Any layer creation options are ignored as it is an append.


texts: Problem Docking Add-In Toolbar in ArcGIS 10 SP4/5

texts:
I have developed an Add-In Toolbar for ArcGIS 10.  It docks just fine on my development computer (which is running 10.0, no service packs according to ESRI's patch finder).  It has been deployed to users who are running 10.0, SP4 and SP5, and no matter what they (or I) do, I cannot get the toolbar to stay docked where they want it to dock.  I do not set the docking behavior in the code so I can't see any way I might be inadvertently affecting it programmatically.

This issue seems identical to the one here: http://support.esri.com/en/bugs/nimbus/TklNMDYwNDI2, which was supposedly fixed in SP3.  I wonder if it broke again in SP4???  Anyone else notice this?  My main tester seems to remember that when she first tested it (on SP3) it docked just fine, but at some point after that it stopped docking - looking at the timeline of her testing, there is a good chance the problem started manifesting after she upgraded to SP4.

I have tried every trick I can find online (all to no avail):

"Save all customizations to the document" is unchecked
"Create new toolbars and menus in the document" is unchecked
I have deleted normal.mxt and had ArcGIS recreate it
I have run AdvancedArcMapSettings.exe and changed monitor dimensions to match the actual dimensions of the monitor

Any thoughts or other suggestions?  Thanks!!


texts: How do I make the scalebar accurate when I am working with a google satellite image in QGIS?

texts:
I start a project in QGIS. I set the CRS to anything (in this case I use EPSG=26918 a NAD83 datum). I also enable on the fly projection. Then I import a vector. In this case a kml, which I believe is a CRS of WGS84. Then I use openlayers and add a google satellite image. Everything looks good and lines up well. I then open a composer window and add a map and a scalebar. When I check to see if the scalebar is accurate I find that it is not.

If I go back to the project and check the measurements with the measurement toolbar I find that the measurement is accurate when the ellipsoid button is checked and inaccurate when the ellipsoid button is not checked. The measurement that I get when the ellipoid button is not checked is the same as I get in the composer window when I use the scalebar to figure out a length.

I guess I have a sense of why this happens. But I want to know how I can get an accurate scalebar in my composer window when I am using a google satellite image?

Thanks,
David


texts: Slope obtained with GDAL has weird lines

texts:
I am working with this DEM obtained using terraincache, it uses mapzen data. The SRC is WGS 84.

I ran this GDAL command to get the slope:

gdaldem slope -s 111120 dem.tif slope.tif (I had to do scaling because of the SRC)

The resulting image looks fine but it's not smooth, it has weird lines over it, like contour lines.

You can download the DEM file here: https://gofile.io/?c=Du7WnM

Does anyone have an idea of how to solve it?


texts:
I have a single band raster with values 1-9 (9 being the best and a cell size of 5m) this raster has been clipped to a feature class that has approximately 80 individual buffered areas of 250m (none of these buffered areas are spatially joined). Within each of the buffered areas of the raster, there are multiple cells with different values.

I need a single point within each buffered area that has the best value, but for the largest concentration of those values.
If there is no value of 9 then a single point for a value of 8 and so on.

So I will end up with 80 points that indicate the highest concentration of the best value available. Just to be clear the point won't necessarily be in the middle of the buffered area, but somewhere within it.

How do I go about modeling this or a python script in ArcGIS Desktop 10.3.1?


texts: Showing infowindows for overlapping polygons?

texts:
I am using Carto.js and I have a map with multiple polygons with the same boundaries.

One layer is purple, and one is orange. I have tried to use transparency to indicate that sometimes there is an overlap - the overlapping polygons are the ones in the middle that are an intermediate colour:

The problem is that each polygon has an infowindow, and when the user clicks on the polygons in the middle, the infowindow only appears for the top layer. So it's not obvious to the user that there are two overlapping polygons there.

How can I show the infowindows for both polygons when the user clicks?

Or is there another, better way to handle this UX problem?


texts: How to add links to popups in Mapbox GL JS

texts:
I have managed to make a map w/ styled pops and such. However, what I am trying to do now is to have it such that the links inside of the popups connect directly to websites. If possible, I would like the title text (h3) of the popups to go directly to the websites when clicked.

I am a bit of a js/css/webmapping newbie.

Code:

    map.on('click', function(e) {
 var features = map.queryRenderedFeatures(e.point, {
layers: ['masterdirectory-qc-v5-dyz490']
 });

 if (!features.length) {
return;

}

var feature = features[0];

 var popup = new mapboxgl.Popup({ offset: [0, -15] })
.setLngLat(feature.geometry.coordinates)
.setHTML('
' + feature.properties.Company + '

' + feature.properties.City + '
' + feature.properties.URL + '
')
.setLngLat(feature.geometry.coordinates)
.addTo(map);

});

My goal is to take all the URL's contained in the column referenced by feature.properties.URL and connect them to the title, retrieved by feature.properties.Company.

My research turned up Display a link in a popup with Leaflet which seems to be a very similair question. I'm guessing my issue is less a problem understanding mapbox gl js and more a lack of knowledge about webdev in general.


texts: Showing huge GIS data on map using Leaflet?

texts:
I have a huge amount of GIS data in a database. I am fetching this data and converting it into GeoJSON before viewing the map. I am fetching all of the data to meet the requirements, which is taking a lot of time. Is there any way to load data for a particular map window which user is viewing, and only loading respective GIS data for the same window using the Leaflet library?

In my case, I am using Google Maps to show parking slots (GIS data) of a country on a map.


texts: Adding LWGEOM to QGIS 3 Processing Toolbox

texts: In 2.x versions of QGIS I was able to add LWGEOM to the QGIS Processing Toolbox and run the Make Valid command. Using OSGEO4W, I've downloaded QGIS 3.2.1 and also downloaded the latest LIBLWGEOM, but I can't see a way to add the LWGEOM tools to the Toolbox. Any ideas what I'm missing?


texts: Exporting ArcMap produced layer file histogram?

texts:
I could simply take a screenshot, or reproduce it in other software, but what I see when I go to layer properties>symbology>classify is perfect for what I need.

The only thing is  it would be nice to have picture with better resolution, so I'm wondering if I could somehow export this graph. I have ArcMap Basic.


texts: Precision of marker being placed on polyline

texts:
I'm attempting to place a marker on a polyline between two lat/lngs at a given percentage. Now, I know I can do this with google maps, but I need to be able to calculate it in backend code and return the values to be used in an agnostic solution.

My issue is that I have the following calculation to work out the waypoint:

function half($startLatLng, $endLatLng, $decimal){
  return {
      'lat': $startLatLng.lat + ($endLatLng.lat - $startLatLng.lat) * $decimal,
      'lng': $startLatLng.lng + ($endLatLng.lng - $startLatLng.lng) * $decimal
  }
}

(you can see a working here: https://jsfiddle.net/77my4685/4/)

My only problem is that the precision of the marker placement is out, often by a little bit, sometimes by a lot. If you zoom into the marker on the example you'll see my issue, or see the image below:

I assumed it was the precision of the lat/lngs but including more decimal places doesn't resolve the issue. Does anyone know why this issue is occurring and how I can go about solving it.


texts: Clip Basemap Layer in QGIS or set a limitation of extens

texts:
I have put an OpenStreetMap layer in QGIS and I don't want to zoom or to pan to all the countries there, but I want to focus just in one place.

Is there any possibility to clip the Basemap, or to add some coordinates to limit the view in Qgis?


texts: How to make feature display after converting Features to Graphics?

texts:
I used the "Convert Features to Graphics" tool on a polygon layer.  It worked correctly with my feature no longer shown but a graphic being displayed instead.

However now, I'd like to display my feature again.  The feature layer is turned on and the symbology is set, just as it was before.  I tried deleting the Graphic group from the Annotation Groups in the Data Frame Properties, and I also tried deleting the graphic, with no success.

How do I enable my layer to show on the map again?  I'm using ArcMap 10.2.

If I add the shapefile again using the Add Data button, it appears in the map.  I could use this workaround but it seems like there should be a way to salvage my existing layer, rather than having to abandon a layer each time you convert its features to graphics.


texts: Qgis 2.18.1 text under a degree angle (Grass 7.0.5)

texts:
I have a points layer with an attribute table with in one column the text (in this case streetnames) and in another column the angle in which the text should be to match the road (which is in another layer).

Is there a function in Qgis or a GrassTool to get my text in the right angle?

A screenshot from the relevant columns in my attribute table can be seen here: http://imgur.com/a/WkWr2


texts: Matlab script - how to access from/run in GRASS GIS?

texts:
I am currently working on shell script in GRASS GIS and Ubuntu. At some point, I must run a calculation (based on an ascii grid) that I have available as Matlab script and .exe (created from Matlab I believe).
For that I must exit GRASS (and boot into Windows), but it would be neat if I could somehow run this from inside GRASS GIS (or at least stay inside Ubuntu). For the non-Matlab and Python-expert...what would be a sensible approach? From what I have understood, converting to Python might be the only real option.


texts: Need to convert feature classes regularly deleting/overwriting existing shapefiles and breaking the script to terminate

texts:
I know a little Python/arcpy, enough to barely get me by. My goal was to, on a regular task schedule, convert feature classes to shapefiles and overwrite existing shapefiles, but I realized that arcpy.overwriteoutput = True did not work. So, as a work around, I added an os remove code to remove all files ending in .shp, .dbf. prj etc. I ran it in ArcCatalog and it was successful in deleting some file types but not all of the file types associated with the shapefile, and although I commanded it to delete .shx, .cbg, and .sbn files, they still remained in the folder. I also tried arcpy.CopyFeatures_management to copy the feature classes to the existing shapefiles in the workspace, but I was unsuccessful in even getting a version of a script to work.

Here is an example of a code that works but runs in an infinite loop.

(I replaced the workspaces, folders, and destinationoutput) Please focus on the infinite for loop.

 import os
 import arcpy
 dir_name = C:\Users\Employee\Documents\ArcGIS
 location = os.listdir(dir_name)
 extensions = (".shp", ".prj", ".dbf", ".xml", ".cpg", ".sbn", ".sbx" ".shx")
 for item in location:
     if item.endswith(extensions):
         os.remove(os.path.join(dir_name, item))
 arcpy.overwriteoutput = True
 arcpy.env.workspace = C:\Users\Employee\Documents\ArcGIS\NSD_exampleschemaCOPY.gdb
 destinationOutput = C:\Users\Employee\Documents\ArcGIS\HansenTest
 datasetList = arcpy.ListDatasets()
 waterfcList = arcpy.ListFeatureClasses("", "", "WaterDistribution")
 for datasets in datasetList:
     for fc in waterfcList:
         if fc in waterfcList:
             arcpy.FeatureClassToShapefile_conversion(fc, destinationOutput)
         elif fc == "MST_Temp_Pipeline_Project":
             break
         else:
             arcpy.AddMessage("Complete")

I also tried this to work but it keeps processing with no end in sight:

import arcpy
arcpy.env.overwriteOutput=True
arcpy.env.workspace = r'C:\Users\Employee\Documents\afgdb.gdb'
outWorkspace = r'C:\Users\Employee\Documents\Employee'
datasetLists = arcpy.ListDatasets()
fcList = arcpy.ListFeatureClasses()
cases = ("watertest", "watertest1", "watertest2", "sewertest", "sewertest1", "sewertest2", "sewertest3")
for datasets in datasetLists:
    for shapefile in fcList:
        for shapefile in cases:
            outFeatureClass = os.path.join(outWorkspace, shapefile.strip(".gdb"))
            arcpy.CopyFeatures_management(shapefile, outFeatureClass)

I for some reason couldn't answer my own question, but the script below works without looping infinitely.

import arcpy
arcpy.env.overwriteOutput=True
arcpy.env.workspace = r'C:\Users\bnamsaly\Documents\New File Geodatabase.gdb'
outWorkspace = r'C:\Users\bnamsaly\Documents\Vicky'
datasetLists = arcpy.ListDatasets()
fcList = arcpy.ListFeatureClasses()
cases = ("watertest", "watertest1", "watertest2", "sewertest", "sewertest1", "sewertest2", "sewertest3")
for shapefile in cases:
    try:
        outFeatureClass = os.path.join(outWorkspace, shapefile.strip(".gdb"))
        arcpy.CopyFeatures_management(shapefile, outFeatureClass)
    except arcpy.ExecuteError:
        tb = sys.exc_info()[2] # execute errors from Python
        tbinfo = traceback.format_tb(tb)[0]
        pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
        arcpy.AddError(pymsg)
        msgs = "ARCGIS ERRORS:\n" + arcpy.GetMessages(2) + "\n"
        arcpy.AddError(msgs)


texts: Cannot access GeoServer on http://localhost:8080/geoserver

texts: I'm getting error 503 not found for localhost:8080/geoserver/web I downloaded GeoServer version 2.13.2 and pointed the JRE to here: C:\Program Files\Java\jre10.0.1


texts: Only able to view 2000 records of Excel table using ArcMap

texts:
I have added a 97-2003 Excel document to ArcMap to compare data. The file has approximately 3200 records. When I initially open the table in Arc it shows all 3200 records, but as soon as I attempt to sort on a field or even select one of the records the table suddenly only has 2000 records. Some have suggested to go to bottom of the table and select the arrows to scroll to the last record. However, when I do this, the table selection goes to record 2000; the table really does only have 2000 records. If I close out and view the attributes again all 3200 reappear. I did seem to remedy the situation by creating an empty personal GDB and importing the table but I feel like I shouldn't actually HAVE to do this to display 3200 records.

Has anyone else encountered this issue?


texts: A original GeoJSON file is also deleted when I edit an imported layer by PyQGIS

texts:
I am trying to import a geojson file of county in U.S. and editing it. Using the following code, I can get a features of the state (in this case "Missouri").

My question is when I edit the layer, features in the original file ("county.json") are also deleted. To avoid this, how should I do?

QgsProject.instance().clear()

county = QgsVectorLayer('[some path]/county.json',
                        'county',
                        'ogr')

QgsProject.instance().addMapLayer(county)

with edit(county):
    request = QgsFeatureRequest().setFilterExpression('"STATE"!=\'%s\'' % 29)
    request.setSubsetOfAttributes([])
    request.setFlags(QgsFeatureRequest.NoGeometry)
    for f in county.getFeatures(request):
        county.deleteFeature(f.id())


texts: Outputting bookmarks in TXT file from ArcMap using ArcPy

texts:
The only output format permitted by the ArcMap bookmark manager (v10.7) is an "ArcGIS place file (*.dat)" which is not entirely human-readable when opened in a text editor.
What is the file structure of an ArcGIS Place File? suggests it may not be possible to access content directly and probably requires ArcObjects to extract bookmark info.
Is there another way to transfer bookmarks from ArcMap so that they are suitable to be read by QGIS?


texts: How to Calculate the Maximum Distance Between Features?

texts:
I'd like to calculate the maximum distance between polyline features.  I have a base polyline with many other polylines either touching or near the base polyline as shown in the picture.  The distance I would like to calculate is from the base feature (red) to the other features (blue).  Because I'm looking for maximum distance, it would be the perpendicular distance from the "peaks" of the blue polylines to the red polyline.

I'm using FME for this process and I can't seem to find anything that calculates maximum distances.  I thought about chopping the blue lines to points using the "Chopper" transformer and calculating the distance to the base feature and then grouping by some ID and finding the maximum within that grouping.  Then, possibly, put that value back into an attribute of the unchopped line.  However, I'm not really sure how to do all that.

Am I on the right track?  Can anyone advise of another method?  Or how I can implement my method?

Thanks.


texts: Python Script for Zonal Statistics As Table results empty Table

texts:
I have a TIFF format file and want to automate some tasks. I have a Python script for Zonal Statistics as Table but it results in empty table output. However, the Zonal Statistics As Table Tool at ArcMap runs successfully! What is the problem with the Script?
import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension('Spatial')

arcpy.env.workspace = "D:/Data/Hydro/"

raster_list = arcpy.ListRasters("*.tif")
lisTable = arcpy.ListTables()

zone = "D:/Data/Catch.shp"

# Loop through each file and perform the zonal statistics processing
for raster in raster_list:
    raster_name = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
    outTable = raster_name + ".dbf"
    arcpy.gp.ZonalStatisticsAsTable(zone, "NAME", raster, outTable, "NODATA", "MEAN_STD")


texts: How to exclude some classes from legend?

texts:
I have a raster file with 8 classifications (all numerical), however the last classification I have coloured with 'no colour' so it does not show up on the map, however I also wish to exclude this class from the legend. Is there a way to remove certain classes from the legend?

I am using ArcGIS for Desktop 10.5


texts: Exporting ground returns (.las) from classified point cloud?

texts:
I have lidar data that is classified as several categories (water, ground etc.).

I need to open this lidar data in some software (it will be great if it would be a open source software) and delete a couple of them to leave only ground category.

Afterwards I have to export this data as *.las format again.


texts: How correct is the "Accuracy" value given by GPS devices?

texts:
I have a Garmin Forerunner 305 and a Google Nexus One Android phone. From the GPS on both these devices I can get a "Accuracy" value that is usually between 6-20 meters. How much can I trust this value? Is it close to correct or may it be very incorrect? I think it must be derived from the NMEA strings received.

I think that both my devices is using standalone-GPS and Assisted GPS from the mobile phone network, if that matters.


texts: Install PostGIS and TIGER Data in Ubuntu 12.04

texts: Could someone write a brief albeit dumbed down idiot's guide to installing postgis and loading national Tiger data on ubuntu? I've tried a few guides, namely, http://wiki.bitnami.com/@api/deki/pages/302/pdf, but I'm not having much luck. I apologize for the open ended nature of this question.


texts: Return variable with GRASS GIS on a Windows machine

texts:
How do I return a variable in GRASS GIS on a Windows machine?

For example I set a list of files to define my region. I already figured out that I can't use the bash command export, I need to use set instead. The command should be okay, but it doesn't return $MAPS.

set MAPS=`g.mlist type=rast sep=, pat="*"`
g.region rast=$MAPS -p

I tried already to check if the variable is set (echo $MAPS), but without success.


texts: In OpenLayers displaying two images at the same time

texts:
Within OpenLayers, I have two images, both jpg with world files, and each 1440 x 900 px.  The two images are geographically located next to each other.  Is there a way to display both images at the same time.  I can make an image layer for each and use the Layer Selector to toggle between them, but I cannot figure out how to display both images together.

From my reading it appears that an image layer is a base layer, and since two base layers cannot be active at the same time, I can only have one or the other active.

thanks


texts: Can the QField Android App integrate with Dropbox (in order to open & edit/view a desktop Project)?

texts:
I have a project that I work on in QGIS (desktop), for which my project file and shapefiles / layers are all saved in my 'Dropbox' folders.

I would like to be able to open, view and preferably edit the above project in QField on my Android phone when out and about.

I have the Project file and shapefiles / layers visible in my Dropbox App, but the QField App seems unable to integrate with / see my Dropbox folders from which to open up the project.

Perhaps this question should be directed at Developers, or self built... (which I can't do) any suggestions?


texts: ArcPy unzip geodatabase error code 80004005

texts:
I got a zipped geodatabase from a colleague and want to unzip in order to load it to my project. Therefore I run the following command in the ArcGIS Pro Notebook:
import arcpy
arcpy.env.workspace = r"E:\Dateien\ArcMap"
arcpy.UncompressFileGeodatabaseData_management(r"Geodat.gdb")

Unfortunately I get the following error message:
ExecuteError:
 file C:\ArcGIS\Geoprocessing\Geoprocessing\Compress\GPCDFFunctions\COM_utils.cpp,
   line 69: Unbekannter Fehler;  (code 80004005)

("Unknown error; (code 80004005)")

The named filepath ("C:\ArcGIS\Geoprocessing\Geoprocessing\Compress\GPCDFFunctions\") does not exist on my filesystem but I installed ArcMap on another drive (E) and I can't find any file "COM_utils.cpp" on my system.
Is my ArcMap installation corrupted or is there another common reason for the error?
System Information: Windows 10, ArcGIS Pro 2.6.0
I don't use Windows/ArcGIS that often. I'm more a Unix/open source user so I have no clue what is going on.


texts:
I am working on a small QGIS Python plugin and want to get the X/Y coordinates from a user click on the canvas. I have a small UI called from button in a toolbar. After filling the UI with information I would like to get the XY coordinates. I am using QGIS 3 with Python 3 and QT5.

I tried to use the steps from
How to get co ordinates of points on mouse click in PyQGIS
but when running it jumps right to the end without calling the display_point function and there is no reaction to mouse clicks.
Are there any changes to the current QGIS version?

My code, which is inside the plugin class and called by a button from the toolbar without UI.

def run2(self):
    # a reference to our map canvas
    canvas = iface.mapCanvas()

    # this QGIS tool emits as QgsPoint after each click on the map canvas
    pointTool = QgsMapToolEmitPoint(canvas)
    # Checkpoint
    print("S 1")
    pointTool.canvasClicked.connect( self.display_point )
    canvas.setMapTool( pointTool )
    # Checkpoint
    print("S 3")

def display_point( pointTool ):

    print('({:.4f}, {:.4f})'.format(pointTool[0], pointTool[1]))
    # Checkpoint
    print("S 2")


texts: Defining a required output field parameter for an existing or a new field

texts:
I need to define a required field parameter, and I'd like it to work with both existing and new fields, i.e. regardless of whether the field already exists.

The problem is, if I enter a new (non-existing field), I get error 000728: Field ... does not exist within table.

This must be due to the parameterDependencies property, but without it I would not get a dropdown allowing to select an existing field. The direction parameter does not seem to matter (Input or Output).

Basically, I'd like to be able to use the drop-down field selection when the field exists, and enter the name directly for a new field. This sounds pretty basic, I must be missing something simple. Could this be done?

Here is a sample *.pyt code to illustrate what I am trying to do.

import arcpy

class Toolbox(object):
    def __init__(self):
        self.label = "[Test] Select existing or new field"

        # List of tool classes associated with this toolbox
        self.tools = [Test_SelectExistingOrNewField]

class Test_SelectExistingOrNewField(object):
    def __init__(self):
        self.label = "Select existing or new field"
        self.description = ""

    def getParameterInfo(self):
        # Feature class where an existing field resides, or a new field will be added to
        p0 = arcpy.Parameter(
                        displayName   = "Input feature class",
                        name          = "fc",
                        datatype      = "DEFeatureClass",
                        parameterType = "Required",
                        direction     = "Input")

        # Name of the output field in the `normals` feature class which will contain the final (conflated) Z values.
        p1 = arcpy.Parameter(
                        displayName   = "Existing or new field name in the input feature class",
                        name          = "fld",
                        datatype      = "Field",
                        parameterType = "Required",
                        direction     = "Output")
        p1.parameterDependencies = [p0.name]
        p1.value = 'AField'

        return [p0, p1]

    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):
        pass


texts: ArcMap - Opening the attribute table gives 'Could not load data from the data source' error

texts:
ArcMap wont display the attribute table of a dataset.

I get an error 'Could not load data from the data source' and it says 'An invalid SQL statement was used'. It then opens a table with the correct number of fields, but no field headers and no data.

It does however display the attribute table if I make a selection and ask it to show only the selection, and I can select every record in the table and it works.

It also displays the values if I use the identify button.

The datasource is a table in a Geodatabase and was created from FME. FME reports no errors in loading the data, and loads all of the records it is given. Using ArcGIS 10.1 SP1.

The table has no special characters in the table name, or the field names.

Edit 10/05/2013

The Geodatabase will open in QGIS. If I export to another features class in a file geodatabase using ArcCatalog it the works. I have tried using the FILEGDB and GEODATABASE_FILE writers in FME and the attribute table will not open in ArcMap. Both work in QGIS 1.8


texts: Export or download feature layer table data from arcgis online using python or Arcpy

texts: How would I export feature edit layers or their table data from arcgis online using Python or Arcpy?  I want to download the attributes of a feature edit layer.


texts: Reading Oracle data with PyQGIS gives "identifier is not long" error

texts:
Do you have any idea how load Oracle spatial data to QGIS with using the PyQGIS? Is it possible? I try:

uri = QgsDataSourceURI()
uri.setConnection("*host*", "1521", "*database*", "*username*", "*password*")
sql = "(SELECT fid, geom FROM target.ms_ft_bh095_v)"
uri.setDataSource("",sql,"GEOM","","fid")
layer = QgsVectorLayer(uri.uri(), "test", "oracle")

The error message is: ORA-00972: identifier is not long. But FID max value is cca 20 000. Is Oracle driver correct with this syntax? Ms_ft_bh095_v is view in this moment.


texts: Enabling Data Interoperabilty Extension in ArcGIS Pro with Esri Evaluation Account?

texts:
I just created an Esri evaluation account to test the ArcGIS Data Interoperability extension in ArcGIS Pro, I assigned the extension in ArcGIS Online to an account. But I can not get access in ArcGIS Pro.

How to enable the Data Interoperability extension in ArcGIS Pro?


texts: Selecting CRS with different code in QGIS?

texts:
In QGIS (3.4), when you select the CRS EPSG:21036 Arc 1960 / UTM zone 36S, you default to the general one for Kenya/Tanzania (code 1122) with
+proj=utm +zone=36 +south +ellps=clrk80 +towgs84=-160,-6,-302,0,0,0,0 +units=m +no_defs .

However, there are three other codes for this CRS EPSG:21036:

Burundi: Code 3998: +proj=utm +zone=36 +south +ellps=clrk80 +towgs84=-153,-5,-292,0,0,0,0 +units=m +no_defs
Kenya onshore: Code 1284: +proj=utm +zone=36 +south +ellps=clrk80 +towgs84=-157,-2,-299,0,0,0,0 +units=m +no_defs
Tanzania onshore: Code 1285: +proj=utm +zone=36 +south +ellps=clrk80 +towgs84=-175,-23,-303,0,0,0,0 +units=m +no_defs

which of course provide different results (source: https://epsg.io/21036)

Is it possible within QGIS to get to choose between those actually different CRS (an option for those CRS to be shown in the CRS Selector window), or is it necessary to create Custom Projections to use them?


texts:
I'm trying to build QGIS from sources.
I run the ctest command at the end, which gives me such results (master branch):

91% tests passed, 45 tests failed out of 490

Total Test time (real) = 1364.33 sec

The following tests FAILED:
      7 - ProcessingQgisAlgorithmsTest (Failed)
      8 - ProcessingGdalAlgorithmsTest (Failed)
      9 - ProcessingGrass7AlgorithmsImageryTest (Failed)
     10 - ProcessingGrass7AlgorithmsRasterTest (Failed)
     11 - ProcessingGrass7AlgorithmsVectorTest (Failed)
     12 - ProcessingOtbAlgorithmsTest (Failed)
     13 - ProcessingSagaAlgorithmsTest (Failed)
     31 - qgis_coordinatereferencesystemtest (Failed)
     39 - qgis_distanceareatest (Failed)
     44 - qgis_expressiontest (Failed)
     55 - qgis_geonodeconnectiontest (Failed)
     75 - qgis_layoutmaptest (Failed)
     82 - qgis_layoutpicturetest (Failed)
     86 - qgis_layouttabletest (Failed)
    109 - qgis_ogrutilstest (Failed)
    127 - qgis_rasterlayertest (Child aborted)
    144 - qgis_taskmanagertest (Child aborted)
    155 - qgis_ziplayertest (Failed)
    192 - qgis_valuerelationwidgetwrapper (Failed)
    202 - qgis_alignrastertest (Failed)
    207 - qgis_wcsprovidertest (Failed)
    217 - qgis_qgisapp (Child aborted)
    225 - qgis_maptoolidentifyaction (Child aborted)
    248 - PyQgsAnnotation (Failed)
    275 - PyQgsDistanceArea (Failed)
    278 - PyQgsEllipsoidUtils (Failed)
    310 - PyQgsImageCache (Failed)
    320 - PyQgsLayoutAtlas (Failed)
    322 - PyQgsLayoutExporter (Failed)
    335 - PyQgsLayoutLegend (Failed)
    374 - PyQgsOGRProviderGpkg (Failed)
    388 - PyQgsPostgresDomain (Failed)
    399 - PyQgsRasterFileWriter (Failed)
    413 - PyQgsScaleWidget (Failed)
    416 - PyQgsSvgCache (Failed)
    418 - PyQgsTabfileProvider (Failed)
    428 - PyQgsSpatialiteProvider (Failed)
    441 - PyQgsSymbol (Failed)
    443 - PyQgsTaskManager (Failed)
    450 - PyQgsVectorFileWriter (Failed)
    456 - PyQgsVectorLayerUtils (Failed)
    458 - PyQgsVirtualLayerProvider (Failed)
    463 - PyQgsWFSProviderGUI (Failed)
    465 - PyQgsLayerDependencies (Failed)
    468 - PyQgsDBManagerSpatialite (Failed)
Errors while running CTest

I wonder where and how to report these or how to get some information, for example whether or not these failures are of concern?


texts: Trying to Load a Contact in Metadata through Python and ArcMap

texts:
Attempting to create a metadata contact and then use that contact to load in to the Resource Citation Contact and the Metadata Citation Contact metadata fields. The metadata library has an element in the form of:

"citation_contact": {
        "path": "dataIdInfo/idCitation/citRespParty",
        "type": "parent_item",
        "elements": contact_elements},

But I can not seem to properly connect this element "citation_contact" to the new contact that I defined earlier in the same script. In that script, I define a new contact as:

contact = metadata.point_of_contact
contact_name = metadata.point_of_contact.contact_name
metadata.point_of_contact.contact_name = "Bob Smith"

How do I load "Bob Smith" into either Citation Contact field?

Here is the element in the library for the contact:

"role": {
        "path": "role/RoleCd",
        "type": "attribute",
        "key": "value",
        "values": [("resource provider", "001"),
                   ("custodian", "002"),
                   ("owner", "003"),
                   ("user", "004"),
                   ("distributer", "005"),
                   ("originator", "006"),
                   ("point of contact", "007"),
                   ("principal investigator", "008"),
                   ("processor", "009"),
                   ("publisher", "010"),
                   ("author", "011")]},

    "contact_name": {
        "path": "rpIndName",
        "type": "string"},


texts: Why my query returns record type instead of original data type in postgresql?

texts:
I am new to PostGIS and PostgreSQL. I am trying to find out a column of geometry type data and want to makeline or want the text format of the data. So, I used the following queries:

select st_makeline(the_geom)
from (
  SELECT edge.start_geom
  FROM public.edge_table edge
  INNER JOIN routing_result route ON edge.start_id = route.node
  UNION
  SELECT edge.end_geom
  FROM public.edge_table edge
  INNER JOIN routing_result route ON edge.end_id = route.node
) as the_geom

and

select st_astext(the_geom)
from (
  SELECT edge.start_geom
  FROM public.edge_table edge
  INNER JOIN routing_result route ON edge.start_id = route.node
  UNION
  SELECT edge.end_geom
  FROM public.edge_table edge
  INNER JOIN routing_result route ON edge.end_id = route.node
) as the_geom

But both the queries are returning following error:

ERROR:  function st_makeline(record) does not exist

and

ERROR:  function st_astext(record) does not exist

And if I run without st_makeline or st_astext, i.e. the following:

SELECT edge.start_geom
FROM public.edge_table edge
INNER JOIN routing_result route ON edge.start_id = route.node
UNION
SELECT edge.end_geom
FROM public.edge_table edge
INNER JOIN routing_result route ON edge.end_id = route.node

Then it is executing perfectly and returning the geometry values as "record" type instead of geometry type. So, I understood that the above query is not returning the geometry type so st_astext or st_makeline is not executing. But why is this happening and how to fix it?


texts: Issue when I define the CRS of raster files on R

texts:
When I tried to define the CRS of my raster layers env :
crs(env) <- CRS("+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ")

I got this error:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition

Knowing that I'm using R 4.0.3 on Mac, and I didn't have this issue on the previous versions. Does anyone know how to handle this issue?


texts: Can I create a shortcut for Paste Features As?

texts: Is there any way to create a shortcut to "Paste Features As", or even better "Paste Features As->Temporary Layer"? It's not in the list for Settings->Keyboard Shortcuts.


texts: Getting sync error from Collector for ArcGIS?

texts:
Errors when syncing edits using ArcGIS collector.

I have two users (the only ones using the map) experiencing sync errors when taking a collector map offline and trying to sync up their edits...

Both users downloaded a fresh version of this map for a weekend collection, they then returned with their devices to wifi spots and tried to sync their data, which resulted in the following error.

This function was previously working with this map but has been problematic for the last month. The map is also embeded here if you would like to take a look http://www.clubhubssl.com/rss/CMS1246/

Any suggestions on ways to fix this?

Should I make another map for collections and delete this problem one?

Can I use the same datasets for the new map?


texts:
QGIS 2.18.1 - I've used data-defined labelling (X,Y,Rotation fields) with earlier versions of QGIS.  I can move and rotate labels OK with 2.18.1 but cannot change the label text.  When I click on the 'Change label' button and select a label on the canvas, the Layer properties screen comes up, but all the fields are greyed out and hence uneditable.

Has anyone else had this problem?


texts: Delete the newly created OpenLayers string line

texts:
Create line strings to display the size of the sides of a polygon when selected by the user.
Here you have code:
var array_coordinates = feature.getGeometry().getCoordinates();

for (var o = 0; o < array_coordinates[0].length - 1; o++) {
    var point = array_coordinates[0][o];

    PointA = array_coordinates[0][o];
    PointB = array_coordinates[0][o + 1];

    var dist = calc_dist(PointA, PointB);
    var formated_dist = format_dist(dist);

    var coordinates = [PointA, PointB];
    var line = new ol.geom.LineString([PointA, PointB]);
    dist = line.getLength();

    var line_stile = new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'black'
        }),
        stroke: new ol.style.Stroke({
            color: 'blue',
            width: 10
        }),
        text: new ol.style.Text({
            font: 'bold 11px "Open Sans", "Arial Unicode MS", "sans-serif"',
            padding: [30, 30, 30, 30],
            placement: 'center',
            fill: new ol.style.Fill({
                color: 'black',
                width: 42
            })
        })
    });

    var layerLines = new ol.layer.Vector({
        source: new ol.source.Vector({
            features: [new ol.Feature({
                geometry: new ol.geom.LineString(coordinates),
                name: 'line'
            })]
        }),
    });

    dist = dist.toFixed(2);
    var data = dist + " m";

    line_stile.getText().setText(data);
    layerLines.setStyle(line_stile);

    map.addLayer(layerLines);

}

But I would like when the user clicks on another polygon to uncheck / delete the previously created strings on the previously selected polygon
Explanation:
By clicking the user selects a terrain on the map when he selects another terrain different from the first one and I would like to be able to delete the lines with the size of the thickness from the first selected polygon.

Here I can't give a solution at all :(


texts:
I am looking for a way to spatially quantify the difference between three contour line shapefiles. More specifically, I have two elevation rasters which I intend to create contour lines from, but before doing so I would like to compare them to a reference contour line and see if there is any need to revise the input rasters. I understand that I could go the other way and use the Topo to Raster feature in ArcMap and then compare using the Raster Calculator, however I believe I will have an easier time (and explaining to others) doing so the way I originally intended.

I have ArcMap 10.1 and Surfer 11. Thank you ahead of time.


texts: Running ArcGIS Server rather than ArcGIS Desktop extensions?

texts: What the advantage of running ArcGIS Server extensions instead of running the same desktop extension?


texts: Create a street line feature with right-of-way widths using a street centreline file and property boundary polygon layer

texts:
I would like to create a street line feature that contains the varying ROW width along it's length. I have street centreline and a polygon layer representing the street ROW area. I've managed a workaround to visualize this by generating perpendicular lines at regular intervals along the centreline features, then clipping the perpendicular lines to the street ROW polygon and symbolizing the length attribute.

What I haven't figured out is how to apply this length attribute to the centreline feature itself so I can play with that data while maintaining the attributes of the original centreline file (street name, addresses, etc.) Since the ROW width would vary along the length of each segment, can it somehow be stored as a z attribute in the shape field?

Anyone have any ideas or other approaches to achieve the same result?


texts: how to use gdal_rasterize? how to know which attribute to burn with ogrinfo?

texts:
I have a vector layer I need to rasterize. I guess this is a basic operation, but I haven't found the information I need yet.

First, I have reprojected the vector shapefile:

ogr2ogr -s_srs "+init=EPSG:3763" -t_srs "+proj=utm +zone=29 +north +west +datum=WGS84 +ellps=WGS84 +units=m +no_defs" CLC06_PT_WGS84.shp CLC06_PT.shp

The ogrinfo command on the new shapefile yields the following:

  INFO: Open of CLC06_PT_WGS84.shp'
       using driverESRI Shapefile' successful.

  Layer name: CLC06_PT_WGS84
  Geometry: Polygon
  Feature Count: 46915
  Extent: (248668.309849, 3768747.607570) - (537329.149160, 4365413.617030)
  Layer SRS WKT:
  PROJCS["WGS_1984_UTM_Zone_29N",
     GEOGCS["GCS_WGS_1984",
         DATUM["WGS_1984",
             SPHEROID["WGS_84",6378137,298.257223563]],
         PRIMEM["Greenwich",0],
         UNIT["Degree",0.017453292519943295]],
     PROJECTION["Transverse_Mercator"],
     PARAMETER["latitude_of_origin",0],
     PARAMETER["central_meridian",-9],
     PARAMETER["scale_factor",0.9996],
     PARAMETER["false_easting",500000],
     PARAMETER["false_northing",0],
     UNIT["Meter",1]]
  FID: Real (11.0)

now, I want to convert it to raster using gdal_rasterize.
The original file contains information on land use (classed as attribute data), and I want that information to be rasterized.

I have tried the -3d option, but I was asked for the feature source, which is not the .dbf nor the .shx files created.

Another try was

gdal_rasterize -a FID -sql "select FID, * from CLC06_PT_WGS84" CLC06_PT_WGS84.shp CLC06_PT_WGS84.tif

but I got the error message

  ERROR 4: `CLC06_PT_WGS84.tif' does not exist in the file system

I am using gdal 1.6

How can I achieve this?
many thanks.


texts: Convert GeoPackage to ESRI shapefile format

texts:
I learned that to convert a DWG file with all its layers from AutoCAD to an ESRI shapefile format in QGIS, one has to first import it as a GeoPackage, and then convert it from a GeoPackage to a shapefile. I already did that. But there's no option to convert the GeoPackages to ESRI shapefiles in QGIS. One can only export the gpkg layers to a QGIS layer definition. I have tried installing and using the Qpackage plugin with QGIS already, but it didn't work. I also tried ogr2ogr for the conversion it didn't still work.

I need to convert the layers from my GeoPackage to the same layers in ESRI shapefile formats.

EDITS: I have attached an Image of the viewport showing the GPKG layers I imported from QGIS. Whenever I right-click on the gpkg layer [ME_AREA 1], and navigate to export, it doesn't bring up the option to save to a shapefile. Instead it brings "save as a layer definition" I really don't know what else to try at this stage.


texts: Problem with loading GeoJSON file in OpenLayers generated by PHP

texts:
I'm trying to read some geo-data from PostGIS database using PHP, creating a GeoJSON file and showing it on a base map using OpenLayers.

The PHP code works just fine. The generated data in GeoJSON format is valid (I checked it with jsonlint). When I load the page, the developer tools shows no errors, but also doesn't load any data on the map.

I have just started learning OpenLayers and don't know much about it but saw some examples for this and tried to do the same.

Am I doing something wrong in the Openlayers code?

The code is:

var map = new OpenLayers.Map({
    div: "map",
    layers: [
        new OpenLayers.Layer.WMS(
            "WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0",
            {layers: "basic"}
        ),
        new OpenLayers.Layer.Vector("GJ", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "gj/DB2GeoJson.php",
                format: new OpenLayers.Format.GeoJSON({
                    extractStyles: true,
                    extractAttributes: true,
                    maxDepth: 2
                })
            })
        })
    ],
    center: new OpenLayers.LonLat(-112.169, 36.099),
    zoom: 11
});

and the GeoJson file generated is:

{
"type": "FeatureCollection",
"features": [{
    "type": "feature",
    "id": "97",
    "geometry": null,
    "properties": {
        "popn_total": ["97"],
        "boroname": ["Staten Island"]
    }
}, {
    "type": "feature",
    "id": "66",
    "geometry": null,
    "properties": {
        "popn_total": ["66"],
        "boroname": ["Staten Island"]
    }
}, {
    "geometry": ["{\"type\":\"MultiPolygon\",\"coordinates\":[[[[578620.717363209,4495974.81786636],[578535.357633272,4495915.80999248],[578524.626146389,4495908.44992811],[578475.151658836,4495981.94016235],[578550.53822566,4496033.27562186],[578571.6880873,4496047.75223892],[578620.717363209,4495974.81786636]]]]}"]
}]

}


texts: Mapbox Tileset turns blank at high zoom in Shiny using Leaflet, but not in MapBox Studio

texts:
I uploaded a large (~8GB) geotiff to mapbox using the CLI.

mapbox upload bweinstein.6thBridge_03182020 6thBridge_03182020.tif

Everything appears fine in mapbox studio. Max Zoom is 22 and you can zoom in and see the tiles. It is the top layer in the style. Here is a fully zoomed screenshot.

In Shiny I connect the CARTO integration URL and make sure maxZoomLevels are set

leaflet(data=selected_centroids,options = leafletOptions(minZoom = 7, maxZoom = 22, maxNativeZoom=22)) %>% addTiles(mapbox_url)

Everything is fine at coarse zoom levels

but then it disappears and turns grey as you zoom.

Reading here: https://docs.mapbox.com/help/troubleshooting/blank-tiles/#style-id

For raster tilesets, the uploaded image resolution sets the minzoom and maxzoom levels. Higher resolution images will result in the tileset rendering at more zoom levels.

Note: regardless of maximum zoom level, data can be overzoomed and visualized to zoom 22.

Is this likely to be a problem in leaflet, in mapbox studio, or in shiny?


texts: ArcGIS 10.1 Upgrade Geodatabase check version?

texts:
I am trying to check the version of a file geodatabase before upgrading to 10.1, I have a list of file geodatabases and some of them are already in 10.1. Can anyone suggest me how to update only those with version less than 10.1?

I am using python script to do this:

arcpy.UpgradeGDB_management(r"C:\Test\Base.gdb", "PREREQUISITE_CHECK", "UPGRADE")


texts: Memory errors when doing Zonal Statistics on a multi polygon shapefile in Python

texts:
I want to calculate the mean NDVI with the EthNdvi.tif file for each zone (there are around 650 admin zones) in the the woreda shapefile. I'm having some issues with running the following lines of code:

from rasterstats import zonal_stats
stats = zonal_stats("Eth_Woreda_2013.shp", "EthNdvi.tif")

Specifically it gives me a memory error. I tried to do zonal statistics in QGIS and this works fine. However when I code this in a python script it results in memory errors.

I've uploaded the two files in a google link: https://drive.google.com/open?id=1go78VkkH8MdAyfpGr4Sg6CvI-WU3WgrC if you want to replicate it.

My laptop has 8gb of ram.
Could anybody explain me why this doesn't work in a python script, but it does in QGIS? Moreover does anyone have a tip how I could still script this in python; should I buy more memory for my laptop or should I do something else?


texts: Geopandas: counting the number of raster pixels within a shapefile polygon

texts:
I have a shapefile that I have loaded as a geopandas dataframe, and it has a geometry column containing polygons and multipolygons.

My end goal is to do some zonal computations on a GeoTiff raster, in particular I want to compute the mean value within each polygon, and also count the number of pixels that contributed to that polygon zonal mean. Here's a link to another GIS-SE question of mine that outlines my attempts to do that so far using GDAL and rasterstats.

I am just wondering if geopandas is up to zonal calculations yet? If I open a raster in GDAL/numpy/rasterio and load the shapefile in a geopandas dataframe, is there a way to compute the number of raster pixels inside each polygon?


texts:
I'm trying to make a vector layer factory function of sorts and I'd like the styles to be as flexible as possible so I've followed various examples and arrived at something which should work but it doesn't display anything on the map. Can anyone suggest what's going wrong please? The code is below and I've also created a fiddle.

var allStyle = {
'Point': new ol.style.Style({
    image: new ol.style.Icon(/** @type {olx.style.IconOptions} */({
        anchor: [0.5, 60],
        anchorXUnits: 'fraction',
        anchorYUnits: 'pixels',
        opacity: 1,
        scale: 0.5,
        src: 'assets/img/icon/pin.png'
    }))
}),
'LineString': new ol.style.Style({
  stroke: new ol.style.Stroke({
    color: '#f00',
    width: 3
  })
}),
'Polygon': new ol.style.Style({
  fill: new ol.style.Fill({
    color: 'rgba(0,255,255,0.5)'
  }),
  stroke: new ol.style.Stroke({
    color: '#0ff',
    width: 1
  })
}),
'MultiPoint': new ol.style.Style({
    image: new ol.style.Icon(/** @type {olx.style.IconOptions} */({
        anchor: [0.5, 60],
        anchorXUnits: 'fraction',
        anchorYUnits: 'pixels',
        opacity: 1,
        scale: 0.5,
        src: 'assets/img/icon/pin.png'
    }))
}),
'MultiLineString': new ol.style.Style({
  stroke: new ol.style.Stroke({
    color: '#0f0',
    width: 3
  })
}),
'MultiPolygon': new ol.style.Style({
  fill: new ol.style.Fill({
    color: 'rgba(0,0,255,0.5)'
  }),
  stroke: new ol.style.Stroke({
    color: '#00f',
    width: 1
  })
})
 };

var allStyleFunction = function(feature, resolution) {
var featureStyleFunction = feature.getStyleFunction();
if (featureStyleFunction) {
    console.log('featurestylefunction');
  return featureStyleFunction.call(feature, resolution);
} else {
    console.log('allstyle');
  return allStyle[feature.getGeometry().getType()];
}
};

var vectorSource = new ol.source.Vector();
var vectorLayer = new ol.layer.Vector({
source: vectorSource,
style: allStyleFunction
});
map.addLayer(vectorLayer);


texts: Summarizing individual data in table using ArcGIS Pro

texts:
I am attempting to summarize a a bunch of individual rows in a table I added into ArcGIS Pro, and I'm wondering if there is any type of tool or work flow I could use to do so?

This is kind of a difficult thing to describe in writing, but I basically have 900 rows and 10 columns. There are 70 counties in total, and I am trying to summarize the data in each individual column for the entire county so there will be only 70 rows instead of 900 (but I want to keep all 10 columns). Each row represents a different amount, so I cannot just dissolve to remove redundant info.

Right now I am manually going through the excel spreadsheet and summarizing the data for each county, but it is taking forever and there has to be another way to do this quickly. This is not public data so I cannot upload it, but I will upload an example with generic amounts.


texts: Converting from netCDF to GeoTiff using gdal_translate

texts:
gdal_translate netcdf_file.nc output.tif

I'm trying to generate a GeoTiff using this command, and it works... kind of. The output tif image is just black and white, and I'd like actual colors on the output image. When I do a "gdalinfo" on the netCDF file, it says there is only 1 band, and that the "ColorInterp" is "Undefined."


texts: How to remove first character from featureclass in Python?

texts:
I have a bunch of feature classes in a geodatabase that start with an underscore and want to remove it, but ONLY the first instance.  See below:
e.g. _COGCC_W_S --> want it to be COGCC_W_S

I have a script written to find and replace but it keeps replace all values instead of only the FIRST instance.  I thought fc.replace(oldString, newString, 1) should work but it keeps removing all underscores.

arcpy.env.workspace = arcpy.GetParameterAsText(0)
findString = arcpy.GetParameterAsText(1)
newString = arcpy.GetParameterAsText(2)

fcList = arcpy.ListFeatureClasses()

for fc in fcList:
    if fc.startswith(findString):
        arcpy.Rename_management(fc, fc.replace(findString, newString, 1))


texts: Setting "R Folder" Path in QGIS Sextante

texts:
I have installed R via the R-Project website. I have QGIS installed via OSGEO4W (currently 1.9.0 Master), and have Sextante plugin installed and working with GRASS and SAGA providers. I am now trying to get R to run, but with no luck. I initially installed R to Program Files, but found mention that perhaps the space in the file path was an issue, so I have re-installed to C:\R\R-3.0.1\

What paths do I need to enter the "R Folder" and "RScripts Folder" configuration area in Sextante? Or am I missing something?


texts: Seeking Nigeria administrative boundaries

texts:
I am looking for Nigeria administrative boundaries (0-3) level
especially the 3 rd level which is called  "senatorial districts".

Where can I find these?


texts: Exporting QGIS polygon shapefiles to Trimble Juno or similar

texts: I am considering buying a Trimble Juno or similar hand-held field computer to assist with vegetation survey.  My main need is to be able to see on it shapefile polygons drawn on my home PC in QGIS.  I am told by a salesman that I would require to buy data collection software in order to be able to do this.  Is this necessary or is there an open source workaround?


texts: Problem mixing the google maps and WMS layer

texts:
I have been trying to use both google maps and WMS layers as an option in the layerswitching panel. But when I switch the base layer from WMS to Google Maps layer and then back to WMS,

the bounds of the map do not come right and the map is not displayed properly.
Even the vector overlay features which display correctly at 1st in the WMS, get crowded into small area on the google maps. I think there is some problem with the re-instantiation of the maps variable during layer switch, due to which the map bounds are not correct.

But still on the whole, please tell me how to solve this problem as it would be great to be able to use any map service in a single HTML page.


texts: How to classify shapefile data into groups?

texts:
I have a map with a shp file containing information on land use. The attribute table has over 10,000 rows.

I need to classify this information into 6 groups to represent 6 colors on my map. I know you can re classify by adding on a new column on the attribute table and put each row into the 6 groups, however this would take years to do 10,00 rows.

Is there a quicker way to do this?

Many thanks all


texts: Double (triple, quadruple) size of polygon using ArcGIS for Desktop?

texts:
I have several polygons, measured in acres.  I'm looking to double the area, either by a buffer, or some other method.  I will need to triple and quadruple in the end, as well.

I have tried to use Field Calculator in a new field, mathematically doubling the size, and using it as a field buffer, yet the resulting sizes are not double, but some other fraction thereof.

My map units are Decimal Degrees.  I also have no problems opening a new project to accomplish this task, as this project is quite large and complicated.


texts: Issue in displaying Folium Map inside a function call using Python

texts:
I am using the following code to display folium map:
import folium

map_hooray = folium.Map(location=[51.5074, 0.1278], zoom_start = 11)
map_hooray

Using the above code the map gets displayed perfectly.
But when I am trying to display the map inside a function call, the map is not getting displayed:
import folium

def map_call(a, b):
    print(a,b) ## prints 51.5074 0.1278

    map_hooray = folium.Map(location=[a, b], zoom_start = 11)
    map_hooray # Calls the map to display but map is not getting displayed

map_call(51.5074, 0.1278)

Can someone please let me know how to resolve this issue.


texts: OpenLayers: how to display the character degrees in the label

texts:
I would like to display temperature in label.
So how can I write DEGREE SIGN? to be display as  °
I was trying write:
& # 176; or  & deg; or \ & deg; and nothing.
Can anybody knows?

I greet.


texts: Writing out comma separated values inside a single cell in a spreadsheet using python

texts:
Please excuse the long subject header. I have gotten pretty good at writing data out to csv files but i seem to have stumbled into this block. Anytime a ',' (comma) is used to join values, it shifts over to the next cell in the spreadsheet. What I want to accomplished is for example:

Inside a single cell, write out cat,dog,mouse

Already checked the csv module in python docs and nothing immediately jumped out at me.

any ideas?


texts: Recommended resources for viewshed generation using Grass GIS 7.4

texts:
I am working on a project where I need to create viewshed areas around towers with a radius of 5000 meters.
The underlying DEM I created using free Lidar data.
Uptill now I did all this on a laptop with a core i7 proc and 8GBs of RAM. Pretty powerful but does not seem to be enough to complete the task within a reasonable timeframe.

Therefore I am considering the setup of a couple of virtual machines on AWS and I need input on what you guys out there would recommend in terms of resources: CPU, amount of RAM etc.

The process is the following:

Merge DEM tiles into one huge TIFF that fully contains the DEM within 5000 m of the tower.
Use this DEM (in TIFF) format in Grass 7.4 (r.viewshed) to generate the viewshed area.
Vectorize the result viewshed raster and delete areas that have a DN value of 0 (not visible)

Even a smaller subset (a sector of 10°-20° from north) took half of the night to produce only the viewshed so I am kind of desperate to speed things up no matter what.


texts: GeoExt (ExtJs + OpenLayers): bug in IE9 when try to change map_panel from container A to container B

texts:
Techonology: GeoExt, which combines ExtJS with OpenLayers

I am trying to dynamically change the app layout by changing a map_panel from container A to container B, and vice verse; basically, for each operation (double click panel header in following codes), the original container A will hide and remove this map_panel from it, then the new container B will add this map_panel object as an item and show up. The following code snippet works fine in Chrome15 and FF8, but in IE9, it has following errors:

SCRIPT5007: Unable to get value of the property 'removeChild': object is null or undefined

OpenLayers.js, line 825 character 412

I am wondering what else I need to add or config in order to make this work in IE as well... any advice are much appreciated......

codes (copy&run):

<!DOCTYPE HTML>
<html>
<head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <script type="text/javascript" src="http://extjs.cachefly.net/ext-3.2.1/adapter/ext/ext-base.js"></script>
    <script type="text/javascript" src="http://extjs.cachefly.net/ext-3.2.1/ext-all.js"></script>
    <link rel="stylesheet" type="text/css" href="http://extjs.cachefly.net/ext-3.2.1/resources/css/ext-all.css" />
    <script src="http://www.openlayers.org/api/2.10/OpenLayers.js"></script>
    <script type="text/javascript" src="http://api.geoext.org/1.0/script/GeoExt.js"></script>
</head>
<script type="text/javascript">
var map, layer;

map = new OpenLayers.Map( ); //define map object
layer = new OpenLayers.Layer.OSM( "Simple OSM Map");
map.addLayer(layer);
map.setCenter(
    new OpenLayers.LonLat(-71.147, 42.472).transform(
        new OpenLayers.Projection("EPSG:4326"),
        map.getProjectionObject()
    ), 4
);

Ext.BLANK_IMAGE_URL = 'img/s.gif';

Ext.onReady(function() {
    Ext.QuickTips.init();

    var mapPanel = { //put map object in a Panel, using GeoExt class 'gx_mappanel'
        xtype : 'gx_mappanel',
        id : 'MAP_PANEL',
        map : map,
        zoom: 6
    };

    //define viewport
    var viewPort = new Ext.Viewport({
      layout: 'border',
        items: [
        { //east container for mapPanel
            title: 'Map',
            id: 'mapEast',
            region: 'east',
            layout: 'fit',
            width: 200,
            border: true,
            frame: true,
            hidden: true,
            items: []
        }, { //center container for mapPanel
            title: 'Map',
            id: 'mapCenter',
            region: 'center',
            frame: true,
            layout: 'fit',
            items: [mapPanel]
        }]
    });

    var mapCenter = Ext.getCmp("mapCenter");
    var mapEast = Ext.getCmp("mapEast");

    mapCenter.header.on("dblclick",function(){ //double click header of mapCenter, map goes to east Container, and centerContainer hides
        mapCenter.hide();
        mapCenter.removeAll();

        mapEast.add(mapPanel);
        mapEast.show();
        viewPort.doLayout();
    });

    mapEast.header.on("dblclick",function(){ //double click header of mapEast, map shows in center Container, east container hides
        mapEast.hide();
        mapEast.removeAll();

        mapCenter.add(mapPanel);
        mapCenter.show();
        viewPort.doLayout();
    });
});
</script>
<body>
</body>
</html>


texts: Exporting two vectors layer as one image

texts:
I am writing Python code to take two vector layers as input add different symbology to it and then export that as a image. The vectors layers have been added to the map canvas but while exporting only the second vector file in the image not both vector file.
  import os

  w_layer = r"path to first vector layer"

  o_layer = r"path to second vector layer"

  image_location = r"export image path"

  wlayer = QgsVectorLayer(w_layer, "W layer", "ogr")

  QgsProject.instance().addMapLayer(wlayer)

  wlayer.loadNamedStyle(r"path for qml file to display a particular symbology")

  olayer = QgsVectorLayer(o_layer, "o layer", "ogr")

  QgsProject.instance().addMapLayer(olayer)

  olayer.loadNamedStyle(r"path to qml file to display a particular symbology")

   vlayer = iface.addMapLayers(wlayer , olayer)

   settings = QgsMapSettings()

   settings.setLayers([vlayer])

   settings.setBackgroundColor(QColor(255, 255, 255))

   settings.setOutputSize(QSize(800, 600))

   settings.setExtent(vlayer.extent())

   render = QgsMapRendererParallelJob(settings)

   def finished():

      img = render.renderedImage()

      img.save(image_location, "png")

   render.finished.connect(finished)

   render.start() `


texts: Postgis AddGeometryColumn for WKT with different types i.e. Polygon and Multipolygon

texts:
I have some polygon data

first;1;MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))
second;2;POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))
third;3;POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))

and points

1,30,10
2,30,10
3,40,20
4,40,30
5,20,20
6,21,21

For these

CREATE TABLE polygons (
  id            CHARACTER VARYING(50),
  value            INT,
  geo                 TEXT
);
CREATE EXTENSION postgis;
CREATE TABLE points (
  id  CHARACTER VARYING(50),
  lat DOUBLE PRECISION,
  lng DOUBLE PRECISION
);

Then points are prepared like

COPY points (id, lat, lng)
FROM '/dataFiles/samples/points_sample.csv'
WITH
DELIMITER AS ','
CSV HEADER;

-- Add point geometry column to table
SELECT AddGeometryColumn('public', 'points', 'geom', 900913, 'POINT', 2);

-- Populate column with point geometries
UPDATE points
SET geom = ST_SetSRID(ST_MakePoint(lng, lat), 900913);

-- Create a spatial index on points
CREATE INDEX idx_points
  ON points USING GIST (geom);

VACUUM ANALYZE points;

for the polygons I use

COPY polygons (sitename, db_range,geo)
FROM '/dataFiles/samples/geometry_sample.csv'
WITH
DELIMITER AS ';'
CSV HEADER;

However, when trying to convert to a geometry column I run into the following problem:

SELECT AddGeometryColumn('public', 'polygons', 'geom', 900913, 'WKT', 2);

will not take WKT but rather only POLYGON or MULTIPOLYGON. How should such a column be processed?  To not have the statement below fail?

UPDATE polygons
SET geomW = ST_SetSRID(ST_GeomFromText(geo), 900913);


texts: PyQGIS: using changeAttributeValues on new feature

texts:
I try to write a Python script with PyQGIS to update automatically a field when a feature is added or when a geometry is changed. It's work well with the geometryChanged listener but not with the featureAdded listener : For a new feature the QgsFeatureId is negative before a commit (-2 for example) and it seems that changeAttribute does not like that :

Layer point: Feature -2 for attribute update not found.

How can I update a field automatically before a commit?
Here is my code :
#### variables ####
layerName = "point"
fieldName = "test"

## create layer object
layer = (QgsMapLayerRegistry.instance().mapLayersByName(layerName))[0]

## function updateField
def updateField(QgsFeatureId):
    field = layer.fieldNameIndex(fieldName)
    layer.dataProvider().changeAttributeValues({QgsFeatureId:{field:"updateTest"}})

## listeners
layer.geometryChanged.connect(updateField)
layer.featureAdded.connect(updateField)

I work with QGIS 2.14.4


texts: import rasterio gives DLL load failed: The specified module could not be found error

texts:
I tried pip install rasterio but it failed building a wheel. So i downloaded the wheel file and tried installing it again as suggested in this post:
https://github.com/mapbox/rasterio/issues/603
It worked and rasterio-1.0a12 was installed successfully but when I am importing it in python it gives me this error:
DLL load failed: The specified module could not be found.

As suggested in this site https://pypi.python.org/pypi/rasterio , i added gdal111.dll and gdal-data file path to windows path but the same error persists. Any suggestions?


texts: How to apply log scale when plotting sf data?

texts:
I want to plot some sf data, similar to this example:

library (sf)

fishnet <- st_make_grid ()
n <- length (fishnet)
vals <- rlnorm (n)
fishnet <- st_sf (fishnet)
fishnet <- cbind (fishnet, vals)

c_ramp <- colorRampPalette (c ("azure", "yellow", "red"))
plot (fishnet, pal  = c_ramp)

I want to use a log scale for the legend and color scheme, but the plot function doesn't seem to support log scale out of the box. Is there any painless way around that?


texts: Spatial query, embed in layer or project

texts:
I have 2 vector layers:

1) Polygons with parcel-boundaries

2) Polygon with project-boundaries

I want to query all the parcels that are within the project-area, using spatial query 'select parcels that intersect project-area'.
Works fine, so far no problem.

The parcel-layer however is a layer that is often subject to changes.
Instead of always renewing the spatial query, I want to embed the query in a layer or project. So that every time when I open the layer/project, I automatticaly get the current parcels within the project-area.

Can this be done in QGIS?
In Mapinfo it was very easy but I'm walking over to Qgis, so it would be great not to lose this ability!


texts: Keep track of current layers in PyQt plugin for QGIS

texts:
I have a QGIS plugin which has a combobox (QComboBox) element with a list of layers currently opened in QGIS. The problem is that the layer list is only updated on initialization of the plugin (I do it with QgsMapLayerRegistry method described here), so every layer added afterwards is not listed until I I force plugin reload.

So the question is - how do I implement automatic refresh of current layer list in my combo box? For example, it is clearly implemented in all built-in modules like "Select by location" etc.

There is a workaround here which uses click event on a combobox, but I think there must be a more elegant way.


texts: I currently have several versions of GRASS GIS installed on PC with Windows 7 (6.4.3, 6.4.4, 7.0.4, and 7.1). I cannot find the i.ortho.rectify module. It appears to be missing. I have searched and browsed through all of the modules. Any ideas about what might be wrong?


texts: OGR and OSM: Get lat, lon of a point as field for SQL-Query

texts:
I am trying to process a .osm file in python with ogr, but I have problems accessing the lat and lon attributes of the point layer - they are always null.

I want to do spatial filtering to find the closest neighbor in a certain distance, e.g. like:

sql_lyr = osmData.ExecuteSQL("SELECT * FROM points WHERE lat BETWEEN ... AND ... AND long BETWEEN ... AND ...")

(Inspired by the answer to this question: Python: how to speed up spatial search (nearest point)? )

When I call DumpReadable on a feature in the point layer, I get the following:

OGRFeature(points):392633
  osm_id (String) = 392633
  name (String) = (null)
  ...
  lon (String) = (null)
  lat (String) = (null)
  all_tags (String) = (null)
  POINT (16.3540923 48.2041024)

So I see that the lat and lon are read in for the POINT Geometry, but I have no idea how to use this information in a SQL query.

Is there any way I can access the lat and lon via the POINT Geometry in a query, or can I do any further configuration to the OSMdriver so that lat and lon are read in as fields?

I modified osmconf.ini to enable all_tags, to report all nodes, I  tried enabling all common attributes and I also tried adding lat and long to the list of keys like this (even though I think this is pointless, because lat long are attributes of the node tag itself and not really keys...):

# keys to report as OGR fields attributes=name,barrier,highway,ref,address,is_in,place,man_made,lon,lat


texts: How to use different columns of a csv file data as different layers for weighted overlay analysis in QGIS?

texts:
I have a CSV file like this:

I would like to do a weight overlay analysis for each column of the CSV and present it as a layer of the map of United States of America.

The final result should look like this map:


texts: Interpolate for no data areas in a raster

texts:
I have a raster of the average chlorophyll-a concentration for 1 year (NASA MODIS LEVEL2 1km res). However i am still missing values for the pixels near the coastal zone (really important for our study).
Is there anyway I can interpolate and fill these fairly accurately?

So far i have just been experimenting with focal statistics in arcgis, but not really getting accurate results....


texts: Add Jquery UI Dialog Box to OpenLayers 3 overlay

texts:
I'm trying to create a custom popup overlay using Openlayers 3 and Jquery UI dialog boxes. (I want use the dialog box as a popup window).

I'm a bit stuck though as my Jquery coding isn't good. I was wondering if anyone has done this already, and if so could they point me in the right direction? The code I have so far is:

 var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
    view: new ol.View({
      center: ol.proj.transform([-0.06, 51.55], 'EPSG:4326', 'EPSG:3857'),
      zoom: 13
    })
  });

 map.on('click', function(evt) {

        var popup = document.getElementById('popupContainer');
        popup.innerHTML= "
Hello
";
        var coordinate = evt.coordinate;
        $(function() {
        $( "#dialogTest" ).dialog({
                   //autoOpen: false,
                  minWidth: 200,
                  minHeight: 200
                  position: {'within':popup}
             }
        );
        });

         var overlay = new ol.Overlay({
            element:popup

            });

          overlay.setPosition(coordinate);

          map.addOverlay(overlay);
          popup.style.display = 'block';
              $( "#dialogTest" ).dialog( "open" );

    });


texts: Getting points with prevalence rates into a Space Time Cube for Emerging Hotspot Analysis?

texts:
I am working in ArcGIS Desktop 10.4 (advanced license), and also have ArcGIS Pro 1.4 available to me for visualizing in 3D.

My goal is to take disease data for a large dataset of points from 10 different time periods, append them together, and input these points into the Create Space Time Cube tool so I can run Emerging Hotspot Analysis afterward. Information on the disease dataset:

Originally 10 different point features, all with identical attribute tables. The table has 3 fields, a unique ID, a disease flag (1 or 0), and a date field. I do this for all ten, and then I use the Append tool to put them together. I end up having 1 big appended point feature with about 2 million points (~ 200,000 in each time period). Several of the points are coincident in the appended feature but have different dates, which is what I need for looking at disease hotspots over time.

When I input this appended point feature into the Create Space Time Cube tool,

I choose the disease I want to analyze over time as the analysis field, and then the space time cube asks me what kind of summary statistic I would like to use on it as you can see in the image. If I specify the Neighborhood Distance of the bin to 1 mile and choose to sum the disease field, it's going to count all the 1's in the disease field. If they're all 1's and 0's...that's not prevalence. It's just a count of the 1's. If I choose mean, that might work...it'd take the average of 1's and 0's per bin, and look at it over time. But I'd like more control of what goes into the Space Time Cube and I'm not entirely sure what the mean is doing to all of the points being binned.

When I am doing a point-level hotspot analysis for 1 point in time (as in not for the Space Time Cube) I run Integrate with an XY Tolerance of 1 mile, and Collect Events on the points first where disease = 1, then clear my query and run Collect Events again on the total count of points. Spatial join the selection (where disease = 1, with Min summary statistic checked) to the total points Collect Event results. take the resulting spatially joined feature, add a new field called Prevalence, and divide the selection where disease = 1 by the total count. I put this spatially joined point feature into the Optimized Hotspot Analysis tool, choose the prevalence rate as the analysis field, and the result of that tool gives me point-level hotspots. I want to do some sort of equivalent of the point-level prevalence rate to input into the space time cube.

What I am looking for help on:
If I were to choose a 1-mile Integrate XY tolerance, and follow the same steps I'd normally do on regular hotspots (as mentioned in the previous paragraph), when going through the Create Space Time Cube if I were to choose a 1-mile bin size, how could I ensure only 1 point (per time period) w/ an associated prevalence rate would go into 1 bin, and thus preserve my rates going into the Space Time Cube so it's as close as possible to my normal Optimized Hotspot Analysis approach?


texts: How to display only the negative contour labels?

texts: How can I hide the positive contour labels covering the green area and only display the negative contour labels over the ocean areas? I am using QGIS 2.18.2. The contour file is a separate layer and when I enable labels I cannot hide the positive contours. Refer to the image below.


texts: Measuring width of polygons, which is wider than long

texts:
Within my project I am trying to quantify and compare different submarine landslides with ArcGIS. Therefore, I created a polygon shapefile and mapped the spatial extent of each landslide.
Now I want to calculate next to other parameters their exact width and length.
For landslides which are longer than wide, the Minimum Bounding Geometry tool works really fine. However, some landslides are wider than long and this tool is totally misleading in such a case.

Does anyone know if there is a method/tool to calculate the real width and length?


texts:
I have a shapefile with points, and a shapefile with polygons.
I want to sum the attribute value of all points that are contained in a polygon, and I'd like to do this for each polygon in my dataset.
I succeeded to do that with Join Attributes by Location but what I want is that when I modify the polygons, the attribute values must be updated.

Can anyone help me?


texts: Switch an R script from "MODIS Reprojection Tool" to "ModisSwathTool"

texts:
This R script works with MRT (MODIS Reprojection Tool), and its performance is unbelievable for land data.

The only problem is when you use Atmospheric data from Modis, like "MOD05_L2".

The problem comes from here:

Such HDF files (Atmospheric data), only can open with "ModisSwathTool", but this R script is written to use "ModisTool" which can not open such HDF files.

Anyone here has experience to solve this problem? For example, change the R code within the R script to support and work with "ModisSwathTool"?

I mean, anyone had such problem before? If yes, what he/she did to solve that?

Any comments and helps are welcome.


texts: Interactively Display Values of a Heatmap.js leaflet.js Heatmap Layer

texts:
In this example of the heatmap.js tool, I see how the values of the heatmap is returned by using heatmapINstance.getValueAt():

  demoWrapper.onmousemove = function(ev) {
  var x = ev.layerX;
  var y = ev.layerY;
  // getValueAt gives us the value for a point p(x/y)
  var value = heatmapInstance.getValueAt({
    x: x,
    y: y
  });
  tooltip.style.display = 'block';
  updateTooltip(x, y, value);
};

I also see that a heatmap instance is created, the "heatmapInstance" that we are using getValueAt() on:

var heatmapInstance = h337.create({
  container: document.querySelector('.heatmap'),
  // onExtremaChange will be called whenever there's a new maximum or minimum
  onExtremaChange: function(data) {
    updateLegend(data);
  }
});

I've created a heatmap layer using heatmap.js' (see picture below) leaflet-heatmap plugin, how could I display the value of my heatmap layer at the cursor location?


texts: Clipping multiple AOIs from a lascatalog using a dataframe

texts:
I'm tryin to clip multiple AOIs from a Lascatalog using a dataframe containing centres and radii of various plots.

I can probably loop through the data frame and use lasclipCircle() to clip each plot and write it. I was wondering if there is a more efficient way to simultaneously do this while retaining the plot ID (from the dataframe) for each file, if possible.


texts: Extracting metadata like age and resolution from XYZ tilesets using PyQGIS

texts:
I'm using PyQGIS to develop a standalone application. I need to get metadata from the current visible tile, such as the age and resolution of the satellite image.
Is this possible?
It doesn't necessarily have to be in PyQGIS, it's fine if I have to use an external API or something.


texts: ECW files not visible in qgis browser (or qbrowser as executed in console)?

texts:
version: 2.4.0-1ubuntu2

I'm following a web tutorial at lynda.com about QGIS.
In one of the first videos I see there's this preview of raster files.

There's the .dem files and the .ecw files among the tutorial example files. I can preview the .dem files as in the video, but unlike the authors of the tutorial I cannot preview the .ecw files.

Additionally, when I try to open a .ecw file with qgis (not qbrowser), I get the following message:

  The layer whateverthepathis/552B.ecw is not a valid layer and can not be added to the map

When I open a .dem file instead, no error message is displayed and the layer is added just fine.

Is there something I need to setup in order to be able to preview those files and add their layers to the map?
Is it a bug of this particular version?

I'm just getting started with qgis and with GIS software in general.


texts: Saving select query results (year wise) from PostgreSQL/PostGIS to text files

texts:
I have two spatially enabled tables 'bar' and 'foo' in my PostgreSQL 9.5 (x64) database. Table 'bar' contains 87 rows and 'foo' contains 16,060 rows. For each address location (foo.gid), following query selects (260 rows returned in 100ms execution time) nearest streets (foo.geom) within 50 meters as well as interpolates a point to that street. Then, query results are exported to the desired path using postgres COPY() command.

Query:

Copy (
      Select
           foo.gid As addr_ID,
           bar.geom As streets,
           St_LineInterpolatePoint
               (ST_LineMerge(bar.geom),
                  St_LineLocatePoint
                     (St_LineMerge(bar.geom),foo.geom)
           ) As interpolated_point
      FROM foo
      Left JOIN bar ON ST_DWithin(foo.geom, bar.geom, 50)
      ORDER BY
      foo.gid, St_Distance(foo.geom, bar.geom)
     )
To '~path/my_file.txt' WITH DELIMITER ',';

Table 'foo' also contains two columns 'start_time' and 'end_time' (data type 'timestamp without time zone') depicting an interval of varying years like this:

Start_time              End_time
2003-06-07 00:00:00     2004-09-30 23:59:59
2004-02-03 00:00:00     2005-03-10 23:59:59

I would like to break my above select query results into years wise so that I would have a text file of each year with unique records (query results).

Expected output:

my_file_2003.txt
my_file_2004.txt
my_file_2005.txt

Can someone guide me to achieve my expected output?

EDIT:

I am using PostgreSQL version 9.5 with PostGIS version 2.2 on Windows 7 Enterprise (x64) based machine. Thus, I intend to go for Windows' specific solution.


texts: Order of ways in OpenStreetMap relation data (relation members)

texts:
I imported to SQLServer the map of Europe taken from Geofabrik (in pbf format). I preserved the order of nodes assigned to ways and the order of ways assigned to relations. The issue is the order of ways I see on OpenStreetMap website for relations differs from the data I have in db (so in pbf file). The issue is the ways are not in the order they can be linked together to create chain of points to draw some shape. OpenStreetMap website shows proper order of ways (from top to bottom), when pbf has random order of these ways. I noticed it when I was trying to create multipolygon object from Polish border data. The original shape is here.

The data in the table below show that the way/line on position 982 has common points with ways on position 546 and 547, so the way 982 should be between 546 and 547:

Does it mean that OpenStreetMap website reorders ways to show them in correct order or the export from Geofabrik destroyed the order of ways belonging to relations?

Recently I had the issue with wrong order of points belonging to ways forming some shape. The issue has been described here.


texts: Does ArcGIS API support Geoserver's WMS and WFS services?

texts: I want to know whether ArcGIS API for Java Script support Geoserver's WMS and WFS standards. I have published my data on Geoserver and i want to query my data by WFS in ArcGIS API to develope my appliction.


texts: QGIS incredibly slow when trying to update large tables?

texts:
My point postgresql table has about 500k rows and about 20 columns

People have about 5 fields that they update, the other fields are static

Selecting the points is fairly instant, however using the modify selected features, changing 2-3 fields and then hitting OK, it freezes QGIS for a few minutes. Many users are experiencing the same thing

Is QGIS just unable to cope with large tables or am I missing something in terms of optimization?

Saving changes onto the layer is fairly quick, less than 10 seconds.


texts: Pyqgis Fixed Buffer Creation in Meters with Haversine

texts:
On Qgis 2.18, I'm trying to write a python code that allows me to create a 1km, 2km, 3km, 4km, 5km buffer around a selected point feature. To avoid having to reproject things, I'm using the reverse haversine formula. My workflow is as follows:

Get lat long of selected point feature.
Use reverse haversine formula to get coordinates of new points 1km away from origin at specified angle intervals.
Join new points to form buffer ring polygon feature.
Repeat for 2km, 3km, 4km, 5km.

from math import *
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

layer = iface.activeLayer()
for feat in layer.selectedFeatures():
    geom = feat.geometry()

start_x = geom.asPoint().x()
start_y = geom.asPoint().y()
start_lon = start_x * pi / 180
start_lat = start_y * pi / 180
sides = 64
radius = 6378137.0  # meters
points = []
distance = [1000, 2000, 3000, 4000, 5000]

for i in range(len(distance)):
    dist = distance[i]
    degrees = 0
    while degrees < 360:
        degrees = degrees + 360/sides

        bearing = degrees * pi / 180

        end_lat = asin((sin(start_lat) * cos(dist / radius))+ (cos(start_lat) * sin(dist / radius) * cos(bearing)))
        end_lon = start_lon + atan2(sin(bearing) * sin(dist / radius) * cos(start_lat), cos(dist / radius) - (sin(start_lat) * sin(end_lat)))

        points.append(QgsPoint(end_lon * 180 / pi, end_lat * 180 / pi))

        geometry = QgsGeometry.fromPolygon([points])

Two places where I'm stuck:

When I add this script on the Qgis processing toolbox (Add script from file) and try to run it, I get an error which says "iface is not defined".
How can I save each iteration of the distance buffer as a shp file?