texts: How to select pixel with smallest number in stack raster using Python script in ArcGIS for Desktop?

texts:
I have 6 tiff file in ArcGIS which are identical in term of geometry (same frame, same number of row and column). The pixel value of each tiff file is integer and range from 1 to 4.

So if I overlay them on top of each other, for each pixel I will have 6 DN number (think like an RGB color composite). Now I need an output raster with pixel value that take the smallest number in the stack. So for 2 pixel that have 123456 and 334662 it will take 1 and 2 as output respectively.

I would like to do it in Python but if it can be done via Raster Calculator or ModelBuilder it would also be good.


texts: I have some data (points and rasters) that I would like to be able to convert to and from several vertical datums, including MLLW in coastal Florida. I would appreciate it if people could provide some insight into accuracy, pitfalls, requirements, etc. of using available utilities such as VDatum or corpscon. Is one more user-friendly than the other? Are there any others that I should consider? Thanks!


texts: Update a newly added column of a table using ST_Intersection as a filter

texts:
I would be really grateful to anyone that could help me with the following issue.
I have ever used PostGIS much before - so could be considered a newbie...

My Situation:

I have TWO PostGIS tables:

roads
country_boundary

The roads table is so big (within the application 500MB) that I wish to display it by country_boundary.

I can do this successfully using the following query:

select ST_ASBINARY(roads.geog::geometry) from roads, country_boundary
where country_boundary.adm0_name='Sudan'
and ST_Intersects(roads.geog::geometry, country_boundary.geog::geometry)

HOWEVER it takes too long to display the data (several minutes) note: I have tried most optimization techniques, spatial indexing, (ST_Simplify/ST_SimplifyPreserveTopology) these either remove data or takes just as long.

I thought I had found a simple solution that worked - I add a new column 'adm0_name' (aka country name) to the roads table and populated it using the ST_Intersects query above.

update roads
set adm0_name = 'Sudan'
where ST_ASBINARY(roads.geog::geometry) in
(
  select ST_ASBINARY(roads.geog::geometry) from roads, country_boundary where country_boundary.adm0_name='Sudan' and ST_Within(roads.geog::geometry, country_boundary.geog::geometry)
);

The resulting display query: works fast, and as expected.. >>problem solved<< ??

select ST_ASBINARY(roads.geog::geometry) from roads where roads.adm0_name='Sudan'

UNFORTUNATELY NO...

The ST_Intersects includes a small number of geometries (polyline(s) that are shared with neighbouring countries. This means polylines (roads) crossing the country boundary, cannot be give a unique country name. Note: I originally used ST_Within but this removed data that was shared between neighbouring countries (as expected with features like roads/rivers that do not respect country boundaries).

..... thanks for bearing with me...

I decided I needed to use ST_Intersection because this would avoid the boundary issue - it would 'clip' the road dataset to each country boundary thus avoiding the polyline boundary issue.

This works - I have an intersection query  that displays the data OK

Select ST_ASBINARY (
ST_Intersection(g1.geom1, g1.geom2)) FROM (Select country_boundary.geog::geometry As geom1, roads.geog::geometry As geom2 from country_boundary, roads
where country_boundary.adm0_name = 'Sudan') AS g1;

MY PROBLEM: I cannot find a way to update the roads.adm0_name column with the country name because ST_Intersection seems to returns a geometry, it does not behave like ST_Intersects, where I could update the column.

My attempt to do this is below ---> 26 hours later I have reached a point of head scratching. Could anyone of you guys help me out?

update roads
set adm0_name = 'Sudan'
where ST_ASBINARY(roads.geog::geometry) in
(
 Select ST_ASBINARY (
ST_Intersection(g1.geom1, g1.geom2)) FROM (Select country_boundary.geog::geometry As geom1, roads.geog::geometry As geom2 from country_boundary, roads where country_boundary.adm0_name = 'Sudan') AS g1;
);

However the above statement fails ---> I need to update the roads.adm0_name with country_boundary.adm0_name where (country_boundary "ST_Intersection" roads).

Thanks again --


texts: Setting layer name in composer legend as "Hidden" using PyQGIS?

texts:
Is there a way to hide a layer name in the composer legend using pyqgis as we do from the legend properties (see picture)?

The only way I've found so far is by changing the layer name in the composer legend but it will also change it in the Layers list and I don't want that.

This is the code I used:

        itemLlegend = composition.getComposerItemById('Legend')
        modelindex=itemLlegend.modelV2().index(0, 0)
        itemLlegend.modelV2().setFlag(QgsLayerTreeModel.AllowNodeRename)
        itemLlegend.modelV2().setData(modelindex,'')
        itemLlegend.updateLegend()

Any idea?


texts: Getting active PostgreSQL connection username when stored through authentication configuration in QGIS 3

texts:
I've been trying to get the username of a PostgreSQL connection in QGIS 3.6, using this Python code

layer = iface.activeLayer()
provider = layer.dataProvider()
username = QgsDataSourceUri(provider.dataSourceUri()).username()
print (username)

like in this question Getting active PostgreSQL connection's username as variable in QGIS?

It does give me an empty string and I think it is because the connection uses an authentication configuration, so the username is not available this way. Is there another way ? I found the QgsAuthManager class but as a newbie I really can't figure out what to do with it.


texts: QGIS - Disable ability to download plugins for unauthorized users?

texts: In QGIS desktop, is there a way to disable the ability for unauthorized users to download plugins? Is it possible to change the settings so that the "Manage and Install Plugins" menu item is not available?


texts: Error when stop editing a feature in ArcPy

texts:
I want to make changes in a feature attribute table  that is an edge source of a Network dataset. In the next step, I need to rebuild the Network to make the changes take efect.

The problem is that, when I open the editor session and make changes in attribute table, then when I want to close the session the script crashes.

Here is a sample of my code:

import arcpy
from arcpy import env
arcpy.CheckOutExtension("Network")

gdb = r"C:\GDB\NetworkGDB.gdb"
env.workspace = gdb
env.overwriteOutput = True

edt = arcpy.da.Editor(gdb)
edt.startEditing(False,False)
WHERE = '"OBJECTID_12" = 2'
direcao = "yes" #This value was set by another process that's indeppend on this one... Inclusive I tried a fixed value like in this example

with arcpy.da.UpdateCursor("Network_nd_roads",["OBJECTID_12","osm_oneway"],     WHERE) as cursor:
    for lin in cursor:
        txtLog = "Change value from " + lin[1] + " to " + direcao
        escreverArquivoLog(logFile,txtLog) #function to write a logfile
        lin[1] = direcao
        cursor.updateRow(lin)
        escreverArquivoLog(logFile,"Change suceed!!!")

escreverArquivoLog(logFile,u"Closing edition...")
edt.stopEditing(True) #-----HERE HAPPENS THE ERROR!!------

The error is:

  File: C:\SCRIPT\Python\OneWayRules.py

  Message: RuntimeError: A network source with the specified name don't exist.

  ArcPy Message: ERROR 999999: Fail to execute function. User don't have permission to execute the operation. [N_1_Desc] An error happens
  when executing (BuildNetwork).

But as you can see, I dont call BuildNetwork function yet, I just make changes in the edige source of my Network.

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

As  blord-castillo says, I tryed to use with block in my edition and the startOperation(),stopOperation() inside this block, but the error occurs in __exit__ Attribute. Folow my actual code:

with edt.startEditing() as ed:
    WHERE = '"OBJECTID_12" = 2'
        ed.startOperation ()
        with arcpy.da.UpdateCursor(roadsLayer,["OBJECTID_12","osm_oneway"], WHERE) as cursor:
            for lin in cursor:
                txtLog = "Aterar o valor de " + lin[1] + " para " + direcao
                escreverArquivoLog(logFile,txtLog)
                lin[1] = "yes"
                cursor.updateRow(lin)
                escreverArquivoLog(logFile,"Suceed!!!")
        ed.stopOperation()

The error stil appearing in this code, on exit statement.


texts: how do i create OGRSpatialReference from raster file's georeference C++

texts:
I want to create .shp file using ogr with georeference based on my raster file, so i create an ogr layer using CreateLayer( "point_out", NULL, wkbPoint, NULL ); i want to set the OGRSpatialReference variable from NULL to to my raster file's georefrence, all i know to get raster file georefrence is using GetProjectionRef() function. how do i convert this to OGRSpatialReference ? or there is some other way to done that ?

(note i'm using c++ programming language)


texts: Adding an "M" coordinate (dimension) to an sf object in R?

texts: Is there a way to add an "M" coordinate (dimension) to an sf object in R? As far as I can tell, currently, only the "Z" coordinate can be added using the "st_zm()" function.


texts: Bringing in KML into Google Earth with different Coordinate System

texts:
I have a shapefile (in red outline bellow) that has the following attribute:

Projected Coordinate System: NAD_1927_StatePlane_Wyoming_East_Central_FIPS_4902
Projection: Transverse_Mercator

But when I export the shapefile to KML/KMZ so I can view it in Google Earth, the polygon does not align to where the shapefile appears in ArcMap. I can tell because in ArcMap the road is inside the polygon and in Google Earth it's clearly outside see images below)

I've read that Google Earth uses WGS 1984 as its projection, so I projected the shapefile to WGS 1984, exported it to KML/KMZ and still the polygon does not fall in the right place (see purple polygon in image 2) - the road is still out side of the polygon.
Any thoughts?


texts: Python script tool with multivalue raster input not honouring raster path?

texts:
I am creating a custom Python script in ArcGIS that takes a rasters as a multivalue input and have run into a problem which I don't think is a duplicate (having scoured SE for many hours)

My problem is that I want my script to look at the name of the raster and mosaic it if appears more than once in the list of inputs. Note that the directory path can be different, but the raster name the same i.e. G:\DELIVERABLES\ABC.gdb\RASTER and G:\DELIVERABLES\123.gdb\RASTER

Where I am failing to understand is why it does not seem to preserve the data source (rasterpath) when I use the Counter tool from the collections module. I know I am probably going wrong in terms of how I am indenting here the script but I have script fatigue.

If I have an inputRasters value of:

G:\DELIVERABLES\ABC.gdb\RASTER;G:\DELIVERABLES\ABC.gdb\RASTER_1;G:\DELIVERABLES\123.gdb\RASTER_2;G:\DELIVERABLES\123.gdb\RASTER

The print results of arcpy.AddMessage(rasterpath) will give me a rasterpath value of 'G:\DELIVERABLES\ABC.gdb\RASTER_1'. It's giving me the data source of the third raster.

The print results of arcpy.AddMessage(rasterlist) below are:['RASTER','RASTER','RASTER_1','RASTER_2']

import arcpy, sys,os
import re
import collections
rasterlist = []
inputRasters = arcpy.GetParameterAsText(0)
for rasterpath in sorted(inputRasters.split(';')):
   Ras = str(os.path.basename(rasterpath))
   rasterlist.append(str(Ras))
   rasterlist = sorted(rasterlist)
arcpy.AddMessage(rasterlist)
counter = collections.Counter(rasterlist)
for key, value in counter.items():
   if value > 1:
      arcpy.AddMessage(rasterpath)

Can anyone tell me where I am going wrong?

I have tried varying the indentations to no avail.

Do I have to use split up the rasterpath using a "/" instead of using os.path.basename?


texts: Symbolizing layer with graduated colors by date/time field?

texts:
I've got a layer of polygons where each has a date associated with it.  I want to color the polygons with graduated colors, but I want the colors to start with the newest polygons and fade to the oldest polygons.  So far I haven't been able to find a way to do this in ArcMap (<= 9.3.1) other than to create a new field and calculate some integer or float value based on the date.  Is it possible to use a date or time field directly?

The data is natively a date-time value
There are far too many records to use any unique-value renderer
Converting the date to a numerical field is an answer for a different question.


texts: QGIS Rasterize - Use highest value in cell as burn in value

texts:
I am trying to create raster layers of presence/absence bridges in a project area to use as inputs for habitat suitability models. I am using 15 m x 15 m cell size.

The presence of bridges is represented by 6 m diameter circular polygons following conversion of point data to polygon through buffering). The polygons are coded 1 = presence of bridge, and 0 = absence.

When I create the raster, not all of the bridge locations show up (value 0 returned where I want 1). I assume this because the value that covers the majority of the raster cell is being burned in, which in may cases would be 0 - absence.

I want all cells that contain a bridge to return a value of 1, regardless of the total area of the cell that is covered by the bridge polygon.

How would I specify that the maximum value from the cell be used as the burn in value or what are other ways I could achieve this goal?

I am using QGIS 3.4 primarially but also have 3.10


texts: Publishing a FeatureService in ArcGIS Server

texts:
I have a registered PostGIS DB in ArcGIS Server Enterprise > Site > DataStore.

How to create a FeatureService in ArcGIS Server from my PostGIS DB to create a FeatureLayer in JS app?
How to obtain the URL that shows this example ?

const layer = new FeatureLayer({
  // URL to the service
  url: "https://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer/"
});


texts:
I would like to import qgis.core and utils from my own copy of python.  I do have PYQT4 installed there.

I'd be happy to add to my PYTHONPATH or sys.path.append in my code, where do the qgis modules live (on a Mac)?


texts:
I need to know how to install QGIS 2.0.1 for Windows 8 64 bit so it can recognize and open ESRI .gdb files.

What do i download and what settings do i change?

I have seen Installing File Geodatabase (*.gdb) support in QGIS?
but the information at this link was posted two years ago, and I am looking for an update.


texts: Write Excel values in to the attribute table of a shapefile in QGIS

texts:
I am using QGIS version 2.18.19. I have written a standalone script which converts a DXF file in to polygon shapefile. Now I want to write data into the attribute table of the polygon shapefile from a Excel file which I used as a input to the script.

My code is as follows:

dxf_to_polygon=name
t=table
dxffile=file
fpoly=output vector
outputs_GRASS7V_IN_DXF_1=processing.runalg('grass7:v.in.dxf', dxffile,'1',True,False,True,False,False,True,'200000.000000,700000.000000,200000.000000,700000.000000',0,None)
outputs_GRASS7V_CLEAN_1=processing.runalg('grass7:v.clean', outputs_GRASS7V_IN_DXF_1['output'],1,0.1,'200000.000000,700000.000000,200000.000000,700000.000000',-1.0,0.0001,None,None)
outputs_QGISPOLYGONIZE_1=processing.runalg('qgis:polygonize', outputs_GRASS7V_CLEAN_1['output'],True,True,None)
outputs_SAGAPOLYGONDISSOLVEBYATTRIBUTE_1=processing.runalg('saga:polygondissolvebyattribute', outputs_QGISPOLYGONIZE_1['OUTPUT'],'area','layer','perimeter',False,None)
outputs_SAGAPOLYGONDISSOLVEALLPOLYGONS_1=processing.runalg('saga:polygondissolveallpolygons', outputs_SAGAPOLYGONDISSOLVEBYATTRIBUTE_1['DISSOLVED'],False,None)
outputs_QGISADDFIELDTOATTRIBUTESTABLE_1=processing.runalg('qgis:addfieldtoattributestable', outputs_SAGAPOLYGONDISSOLVEALLPOLYGONS_1['DISSOLVED'],'Req_No',2,30.0,0.0,None)
outputs_QGISADDFIELDTOATTRIBUTESTABLE_2=processing.runalg('qgis:addfieldtoattributestable', outputs_QGISADDFIELDTOATTRIBUTESTABLE_1['OUTPUT_LAYER'],'Surveyor_Name',2,25.0,0.0,None)
outputs_QGISADDFIELDTOATTRIBUTESTABLE_3=processing.runalg('qgis:addfieldtoattributestable', outputs_QGISADDFIELDTOATTRIBUTESTABLE_2['OUTPUT_LAYER'],'Survey_Date',2,20.0,0.0,None)
outputs_QGISADDFIELDTOATTRIBUTESTABLE_4=processing.runalg('qgis:addfieldtoattributestable', outputs_QGISADDFIELDTOATTRIBUTESTABLE_3['OUTPUT_LAYER'],'Plan_Name',2,20.0,0.0,fpoly)

How can I solve this?


texts: QGIS Adding single arrow to polyline with dynamic position

texts:
In QGIS, is it possible to add a single arrow (indicating direction) to a polyline  with a dynamic position, so that it is always in view as long as the polygon is in view?

Adding "Arrow" to a line either produces arrows at each vertex (checked repeat arrow on each segment) or does not follow the curvature of the polyline (unchecked repeat arrow on each segment)


texts:
I'm digitizing some old topographic maps (scanned images), and these maps are in grayscale.

When in editing mode, to draw the contour lines, I would like to be able to change the color and size of the cursor, because the default one is too small, and the black color makes it disappear.

Something like in Blender 3D, where we can define, in pixels, the cursor size. Or in AutoCAD, where we can choose any color to the cursor.

I'm using QGIS 2.10.


texts:
I extracted vertices using 'Extract vertices' in Processing Tool. But so many unnecessary vertices extracted, so I want to remove them. What I want to leave is at the beginning and the end of the line. How can I remove those vertices in the middle of the line? I would like to know if there is a way to erase it all at once because there are so many things like that.
Vertices are in the same layer.

I want to remove all the blue vertices like the picture below.


texts: Checking field type in ModelBuilder using if else?

texts:
I am creating a model in gis ModelBuilder.
I want to check the type of a particular field and

if it is long integer than the tools should be different and
else (i.e. not a long integer) my tools are different.

Is there a way to check the field or any already available python script which could help me do this?


texts: Use QGIS to identify and count pixel colour in a photo

texts:
I am trying to use QGIS to analyse the colour component of a photo. I have imported the photo as a raster layer, created a polygon of the area in question and have tried to use zonal statistics but this doesn't appear to work. When I use the zonal statistics function, I fill in all the relevant boxes and click ok. Then the message window closes. That's it. When I have tried using the function on maps I already have, after clicking ok the progress bar came up and after it finished I could find the resulting stats. When I tried it on my photo, absolutely nothing happens. The window just closes, there are no stats anywhere.

Ideally, the stats I would like to produce are the number of pixels in the polygon, and the number of pixels of each colour. Or, to assign the colours into similar groups say 0-25, 25-50, 50-75 etc up to 255. Of course, this is made more challenging by the fact that a photo used as a raster layer is made up of 3 colour bands and, as yet, I haven't managed to separate the bands.


texts: ERROR: could not load library "C:/Program Files/PostgreSQL/11/lib/rtpostgis-2.5.dll": unknown error 126?

texts:
I have got this message when I wanted to create a new extension.

How can I solve this problem? I have Win10 64 bit and PostgreSQL 11.


texts:
I've made an arcpy script which worked last week, but I placed it somewhere different on my computer. I've made sure all the source files etc. are connected. However now I keep getting an error message when it comes to the erase function, it comes up with the following error:

Traceback (most recent call last):
  File "G:\zn\NM\GGB\Applicaties-Admin\Kerngis\ExpImp\Base data\Scripts\List mxd.py", line 22, in <module>
    arcpy.Erase_analysis (lyr, Export, Import)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\analysis.py", line 176, in Erase
    raise e
ExecuteError: ERROR 999999: Error executing function.
The table was not found.
The table was not found. [Kabel_en_leiding_lijnen_Erased]
The table was not found.
The table was not found. [Kabel_en_leiding_lijnen_Erased]
The operation was attempted on an empty geometry.
Failed to execute (Erase).

Failed to execute (Import)

My script is the following:

import arcpy
import arcpy.mapping

Project_number = arcpy.GetParameterAsText(0)

arcpy.env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument("G:\\zn\\NM\\GGB\\Applicaties-Admin\\Kerngis\\ExpImp\\{}\\{}.mxd".format(Project_number, Project_number))  # Uses your currently open MXD
for df in arcpy.mapping.ListDataFrames(mxd, ''): # Loop through dataframes
    for lyr in arcpy.mapping.ListLayers(mxd, '', df): # Loop through layers

        # Any tools you want to run on each layer go here

        Export = "G:\\zn\\NM\\GGB\\Applicaties-Admin\\Kerngis\\ExpImp\\{}\\Export_{}.gdb\\Mask".format(Project_number, Project_number)
        Import = "G:\\zn\\NM\GGB\\Applicaties-Admin\\Kerngis\\ExpImp\\{}\\Import_{}.gdb\\{}_Erased"(Project_number, Project_number, lyr)
        NewFile = "G:\\zn\\NM\GGB\\Applicaties-Admin\\Kerngis\\ExpImp\\{}\\Import_{}.gdb\\{}_MultipartToSing".format(Project_number, Project_number, lyr)
        Merge = "G:\\zn\\NM\GGB\\Applicaties-Admin\\Kerngis\\ExpImp\\{}\\Import_{}.gdb\\{}_Merged".format(Project_number, Project_number, lyr)

        arcpy.SelectLayerByLocation_management(lyr, "INTERSECT", Export, "", "NEW_SELECTION", "NOT_INVERT")
        arcpy.Erase_analysis (lyr, Export, Import)
        arcpy.DeleteFeatures_management (lyr)
        arcpy.Append_management(Import, lyr, "NO_TEST")
        arcpy.Append_management(NewFile, lyr, "NO_TEST")

I have also tried running these same steps in modelbuilder, to see if it was the script or if something else was going on. The same error occurs, the image below will display my model.

I cannot figure out for the life of me why it is giving me errors. The feature class being used is a line feature class, not sure if that makes a difference.

I am using ArcGIS 10.3.1


texts: How to deal with minus values on canopy height model?

texts: I created a canopy height model (from lidar data) using lidR package and wanted to check it on Arcmap. I see there are minus values on canopy height models like -0.246219, -0.0302621. How can I deal with minus values in my canopy height model?


texts: Why does CON statement give ERROR 010240: Could not save raster dataset to (value) with output format GRID?

texts:
I am receiving the following error when trying to run the following code:

import arcpy

from arcpy import env

from arcpy.sa import *

arcpy.CheckOutExtension = "Spatial"

env.workspace = r"C:\final_mosaic_files\for_python"

Z3 = Raster("mos_zb3")

Z4 = Raster("mos_zb4")

PC2 = Raster("pca_finb2")

outCon = Con(Z4 < -0.9274,2, Con(Z3 <-0.6223,2, Con(PC2 > 1.13,2, Con(Z3 <-0.05342,2,1))))

outCon.save(r"C:\classifications\actual_real_deal\CART")

The error is:

  RuntimeError: ERROR 010240: Could not save raster dataset to
  C:\Users\User\Documents\Mo with output format GRID.

I have tried changing the destination folder, and the file does not currently exist in the directory, and the values of the raster are within GRID parameters for the type of format. The input rasters are mosaics of z-score values of 12 landsat images delineating coastal marsh in the northeast. The cell size is 30X30m, using WGS1984 coordinate system, no projection system. I have tried running the script with .tif file format instead of GRID, which does not work, throwing another error (file Z3 does not exist or is not supported). I have also tried running the script with other, smaller GRID files, which throws the same error 010240.


texts: Implementing ArcSDE cost effectively using PostgreSQL?

texts:
I wish to implement an ArcSDE environment without breaking the bank. I would be wanting SQL databases.

If I use PostgreSQL what are the limitations and or the requirements to make it work?

We have 15 concurrent ArcGIS for Desktop Advanced licenses with up to 20 different users.


texts: arcpy.Append_management ERROR 001156 - overshoot field length issue

texts:
I append features from one class to other (using arcpy.Append_management including FieldMappings). The issue is I overshoot the length of target field. I do wonder why Append_management does not deal with such issue, anyway...

For now I can increase the length of target field which fixed the issue. But I still need somehow check the length of input string. And do sub-string if needed so the tool won't crash whenever there too long string in input data. Any simple way to achieve this?


texts: How to make QGIS processing algorithms use default parameter?

texts:
My Question is related to: How to make Processing algorithms use default parameter values in QGIS 2.14?

I would like to use gdalogr:translate and only set the values of "INPUT","OUTSIZE","EXPAND" and "OUTPUT".

However, when I run:

processing.runalg('gdalogr:translate', {"INPUT":rlayer,"OUTSIZE":j,"EXPAND":2,"OUTPUT":output})

I get the following error (even tough PROJWIN is not a mandatory parameter):

Error: Missing parameter value for parameter PROJWIN.

I would like to let this parameter as default.
When I change the code and include one of the following snippets:

"PROJWIN": ''    OR   "PROJWIN": None

The following error occurs:

Error: Wrong parameter value  for parameter PROJWIN.

Alternatively, I tried setting the paramter PROJWIN to the extent of my current raster layer, without success:

fileInfo = QFileInfo(f)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(f, baseName)
e = rlayer.extent()
processing.runalg('gdalogr:translate', {"INPUT":rlayer,"OUTSIZE":j,"EXPAND":2,"PROJWIN":e,"OUTPUT":output})

Does anyone know how this could work?


texts: Can X-Ray for ArcCatalog handle attribute indexes?

texts:
I've been using the excellent ArcCatalog X-Ray for working with designing a geodatabase.

Is there a way to create attribute indexes in X-Ray?  The geodatabase XML schema does contain an  element, but I can't find any mechanism or guidance if X-Ray can handle indexes.


texts: Algorithm for finding irrregular polygon centroid (label point)

texts:
I need to find a centroid (or label point) for irregularly shaped polygons in Google Maps. I'm showing InfoWindows for parcels and need a place to anchor the the InfoWindow that is guaranteed to be on the surface. See images below.

In reality I don't need anything Google Maps specific, just looking for an idea of how to automatically find this point.

My first idea was to find the "false" centroid by taking the average lat and lngs and the randomly placing points out from there until I find one that intersects the polygon. I already have the point-in-polygon code. This just seems awfully "hacky" to me.

I should note that I don't have access to any of the server side code outputting the geometry so I can't do anything like ST_PointOnSurface(the_geom).


texts:
I have two raster layers with a count and value in the attribute table but when I try to reclassify with count and equal interval quartiles the classification reverts back to "value" in the symbology. I need to be able to get the symbology  to be based on count and not value. This is important because I need to create a new layer such that only the highest quartiles from the 2 layers are shown.

How can this be accomplished?

Regarding previous operations I needed to specify certain land cover features such that water was removed so I used a reclassify and then ascribed nodata to water...

then I used a extract by mask so the DEM layer also removed water.

Then I performed 2 different hillshades with different azimuth and altitudes and these layers are what need to be classified correctly so they can be added.


texts: Can ArcGIS ModelBuilder be exported to Python and run in QGIS?

texts:
Is it possible to construct a geoalgorithm on the ArcMap ModelBuilder, export this entire algorithm to python syntax (without having ANY knowledge of python) and then import this algorithm to QGIS and run it successfully?

I am more familiar with ArcGIS Desktop and unfortunately the QGIS Processing framework presents enormous difficulties when using GRASS commands to process raster files. I for one am not familiar with GRASS, and the QGIS processing framework is unable to take care of the GRASS part for me when it comes to dealing with rasters.


texts: Replacing white or black by alpha channel in raster in QGIS 2.8

texts: I'd like to replace black or white in black and white raster in QGIS 2.8. I perceived that in some other raster where "nan" value automatically occurred and only white or black is being displayed.


texts: Is there a way to get a good human readable error from a COM error in arcEngine?

texts:
I'm trying to perform the following operation:

 IPropertySet addressProperties = reverseGeocoding.ReverseGeocode(point, false);

And I get this error: Exception from HRESULT: 0x80040217

Is there some function to get a human readable form of this error?


texts: Calculate total area of polygons in shapefile using GDAL?

texts:
I have a shapefile in British National Grid projection:

Geometry: 3D Polygon
Feature Count: 5378
Extent: (9247.520209, 14785.170099) - (638149.173223, 1217788.569952)
Layer SRS WKT:
PROJCS["British_National_Grid",
    GEOGCS["GCS_airy",
        DATUM["OSGB_1936",
            SPHEROID["Airy_1830",6377563.396,299.3249646]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["Meter",1]]
cat: Integer (9.0)

Can I use GDAL/OGR to get the total area of all the polygons in the shapefile, in hectares?

I'm wondering if this is possible with -sql, something like:

ogrinfo -sql "SELECT SUM(ST_Area(geom::geography)) FROM mytable" myshapefile.shp

But trying that I get ERROR 1: Undefined function 'ST_Area' used..

I guess I could import the Shapefile into QGIS, add an area attribute to each polygon, and then sum it, but I would much rather use a command line tool if possible.


texts: Error using grass v.clean: " 'utf-8' codec can't decode byte... "

texts:
I try running the v.clean tool using QGIS 3.4. As input I use a vector line layer (I checked the validity and everything is fine), but running the algorithm I get the following error:

Traceback (most recent call last):
File "D:/Programme/OSGeo4W/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 415, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "D:/Programme/OSGeo4W/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 372, in executeGrass
for line in iter(proc.stdout.readline, ''):
File "D:\Programme\OSGeo4W\apps\Python37\lib\codecs.py", line 322, in decode
(result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xf6 in position 81: invalid start byte

Execution failed after 0.50 seconds

Since it is says it is a decoding error I checked the encoding and it is encoded in UTF-8.

Does anybody know how to solve the error?


texts: Ways to speed up ST_DWithin between millions of segments and thousands of polygons?

texts:
I have 3 PostgreSQL tables with the following data:

tab1 => It contains about 5.7 millions of segments (straight lines)

tab2 => It contains 6673 polygons

tab3 => It contains 19620 polygons

Projected Coordinate System (British National Grid)

All 3 tables have a spatial index

I am trying to create a new table with all segments from tab1 that are less than 250 metres away from either tab2 or tab3. In order to do so I am running the following query:

CREATE TABLE public.final_def_sgm AS
SELECT tab1.segment_id, tab1.gb_def_id, tab1.def_id, tab1.geom
FROM public.flood_def_sgm_50m tab1, public.sop_polyg tab2, public.coastal_1000_polyg tab3
WHERE ST_DWithin(tab1.geom, tab2.geom, 250) OR ST_DWithin(tab1.geom, tab3.geom, 250);

I launched this query about 48 hours ago and it is still running. The problem is that I have no idea whether it will be running for two more days or for twenty more days (who knows).

I previously tried to do the same thing by creating 250m buffers around the polygons and tried to select those segments that intersected with the polygons (all in postgres). Very slow too.

Is it normal that it´s taking days to run the query?
Is there any better way to do what I am trying to do?

"explain" throws the following information:

Gather  (cost=3884.82..447434055344.37 rows=824566434 width=60)

  Workers Planned:

  ->  Nested Loop  (cost=2884.82..447351597700.97 rows=485039079 width=60)

        ->  Nested Loop  (cost=0.00..6414604.08 rows=127287635 width=6444)

              ->  Parallel Seq Scan on sop_polyg tab2  (cost=0.00..1864.88 rows=6488 width=1435)

              ->  Seq Scan on coastal_1000_polyg tab3  (cost=0.00..792.20 rows=19620 width=5009)

        ->  Bitmap Heap Scan on flood_def_sgm_50m tab1  (cost=2884.82..3514.43 rows=1 width=60)

              Recheck Cond: ((geom && st_expand(tab2.geom, '250'::double precision)) OR (geom && st_expand(tab3.geom, '250'::double precision)))

              Filter: (((geom && st_expand(tab2.geom, '250'::double precision)) AND (tab2.geom && st_expand(geom, '250'::double precision)) AND _st_dwithin(geom, tab2.geom, '250'::double precision)) OR ((geom && st_expand(tab3.geom, '250'::double precision)) AND (tab3.geom && st_expand(geom, '250'::double precision)) AND _st_dwithin(geom, tab3.geom, '250'::double precision)))

              ->  BitmapOr  (cost=2884.82..2884.82 rows=1143 width=0)

                    ->  Bitmap Index Scan on flood_def_sgm_50m_geom_idx  (cost=0.00..14.70 rows=572 width=0)

                          Index Cond: (geom && st_expand(tab2.geom, '250'::double precision))

                    ->  Bitmap Index Scan on flood_def_sgm_50m_geom_idx  (cost=0.00..12.19 rows=572 width=0)

                          Index Cond: (geom && st_expand(tab3.geom, '250'::double precision))

I have also tried to obtain the same information using a JOIN in the following way:

CREATE TABLE final_def_sgm AS
SELECT t1.segment_id , t2.coastal_1000_id , t3.sop_id
FROM public.flood_def_sgm_50m t1
LEFT JOIN public.coastal_1000_polyg t2
ON ST_DWithin(t1.geom, t2.geom, 250)
LEFT JOIN public.sop_polyg t3
ON ST_DWithin(t1.geom, t3.geom, 250)

Same thing. Query has been running for more than 48 hours now.


texts: Calculating distance between points from table of UIDs [ArcMap]

texts:
I have a .dbf file containing pairs of point feature Unique Identifiers (UIDs), and need to calculate the straight-line distance between each of these pairs.

The XY coordinates of the individual points are part of a separate feature class (also containing the UIDs).

Is there any way that I could calculate the distance between these points without having to create a new table of XY coordinate values? (There are 3085 pairs altogether, so this would be impractical)

EDIT:

.dbf example:
UID0        UID1
40598473    83894847
40638728    59068967
38989786    49872793

Point feature attribute table example:
UID        X        Y
38989768    425667    346273
40598473    425668    346277
40638728    425598    346256


texts: Difference between name and toponym name?

texts:
When I query geoNames for its data, it typically returns the place name in two ways:

1) name and 2) toponymName

What is the difference between them?


texts: Difference between intersect {raster} and gIntersection {rgeos} in polygons - R?

texts:
I am interested in knowing the difference between tools: intersect {raster} and gIntersection {rgeos}** tools in R. I have several overlapping polygons and I am interested in their overlap, but mostly in spatially merging their attribute tables. From my resulting shp splitted in multiple polygons, I want to calculate area of each of polygons using gArea {rgeos}.

Are there any pitfalls I should be aware of? For both intersection processes and for subsequent areas calculation?

For the moment I found intersect{raster} more satisfactory for my purpose as it keep "SpatialPolygonsDataFrame" (SPDF) and even merge polygons data.frames. However many online helps propose instead gIntersection {rgeos} which produces just "SpatialPolygons" not even "SPDF". There is not really difference in computation time.


texts: Split polygon into parts adding attributes based on underlying polygon in QGIS

texts:
I am trying to split one polygon in several parts adding varying attributes depending on the overlaying other polygons.

To illustrate what I am talking about: the yellow ring is a buffer zone around a marine protected area. I'd like to have this polygon split in different parts (but keep them in the same shapefile) where it overlaps with the light green area (the High Seas) and add an attribute for these areas to the buffer polygon indicating they are located in the High Seas (vs. exclusice economic zones, shown in blue).

I have a large number of shapefiles to do this with, if there is a plugin or batch process I'd be very happy to know.


texts: Seeking data for Mapping Pubs in London?

texts: I want to make a map showing pubs in London but don't know how to go about getting the data. I have looked on digimap and downloaded a few things but none seem to be right. any ideas?


texts: Global Shapefile with USA Administrative Only

texts:
I'm trying to build a simple web application to display service coverage information worldwide.  I have been using Datamaps/D3 to build the application.  I downloaded the Natural Earth shapefiles for countries and for states/provinces. I am using GDAL and Topojson tools to get JSON data for the app.

The problem is that I would like to get a global map of countries with state/province (aka administrative) borders for the USA ONLY.

I have been trolling around the web for a few hours with no luck.  I have downloaded QGIS and working through learning how to do this. But, I figured I would post something here as well.


texts: GEE: Indexing a List of Lists inside a property of a FeatureCollection

texts:
I'm running the LandTrendr algorithm in GEE over multiple points of interest, which leaves me with a FeatureCollection having a list of lists in one the properties for every feature.
Out if this list of lists called LandTrendr I need to know whether the second-to-last element of list #3, is equal to 1 (See orange arrow in the picture). How can I index/query this LandTrendr property to get the desired element?
I have tried many strategies to reach this index, but I don't get anyfurther than isolating the FeatureCollection to a single Feature having only the LandTrendr property left (using .first().select(['LandTrendr'])). At this point it refuses to convert into an array or a list using .toList() or .toArray().

I actually wanted to share the FeatureCollection as an Asset, but somehow GEE didn't let me, so I'll be sharing the whole code below.
Code walk-through: After defining the parameters, building the input time-series, and running the algorithm, I reduce the input back to the original points of interest, as defined by this question. Now, the second-to-last element of list no.3 in the LandTrendr property is the desired index to be retrieved.
// #########################
// # Define inputs and parameters
// #########################

//points of interest
var myPoints = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(-60,-20), {'id': 0}),
  ee.Feature(ee.Geometry.Point(140,-20), {'id': 1}),
  ee.Feature(ee.Geometry.Point(-70,10), {'id': 2}),
  ee.Feature(ee.Geometry.Point(80,40), {'id': 3}),
  ee.Feature(ee.Geometry.Point(10,50), {'id': 4}),
  ee.Feature(ee.Geometry.Point(20, -10), {'id': 5})
]);

// load the LandTrendr.js module
var ltgee = require('users/emaprlab/public:Modules/LandTrendr.js');

// define collection parameters
var startYear = 2005;
var endYear = 2016;
var startDay = '06-01'; // MM-DD start summer North Hem
var endDay = '08-31'; // MM-DD end summer North Hem
var index = 'NDVI';
var ftvList = [];
var maskThese = ['cloud', 'shadow', 'snow'];

// define landtrendr parameters
var runParams = {
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: true,
  recoveryThreshold:      0.25,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};

// #################
// # Prepare Input Time-Series
// #################

// Make annual surface reflection collection
var annualSRcol = ltgee.buildSRcollection(startYear, endYear, startDay, endDay, myPoints, maskThese);

// Make annual LandTrendr collection with 1-observation per year
var annualLTcol = ltgee.buildLTcollection(annualSRcol, index, ftvList);

// add preprocessed time-series to the parameters for the algorithm
runParams.timeSeries = annualLTcol;

// ################
// # run LandTrendr
// ################
var LTresult = ee.Algorithms.TemporalSegmentation.LandTrendr(runParams);

// Reduce the LandTrendr result by the aoi (myPoints). Use ee.Reducer.first() to select
// the pixel that intersects the point.
var ltPoints = LTresult.reduceRegions({
  collection: myPoints,
  reducer: ee.Reducer.first(),
  scale: 30
});

print(ltPoints);


texts:
I am trying to add dynamic text to my map layout that displays the current sheet number out of the total number of sheets (e.g. "Page 1 of 15").  I'd like to set it up so that the label changes when I move from one bookmark (Sheet 1) to another (Sheet 3).

I'm aware of the dynamic text feature that uses data driven pages, but not sure that's the correct method to use... it seems like that is geared more towards working with specific shapefile layers within my MXD file.  Is there another way to accomplish what I'm trying to do?


texts: Printing or Exporting Table on ArcGIS Online?

texts:
I have a table that is filtered on my ArcGIS Online map. The data for this map is from a feature layer that I created using ArcGIS Online by importing the features from a shapefile or .csv.

How would I print or export that filtered table?

I thought it would be a simple task, but I cannot locate any information on the subject.


texts: ArcGIS: Scalable labels in x direction only

texts:
I'm trying to label a road layer that is named in 2 different ways:

XYZ
String between 5 to 20 characters

Is there a way using ArcGIS 9.3 to produce a scalable layer label only in the x direction? Specifically I'm looking for a rectangle with semi circles on either end.


texts:
I work for a company that owns agricultural land.
Our new project is consolidating the small plots we own by swapping with other owners.
We have a village were we managed to find all big owners and to map all our plots. All the owners have plots dispersed all over the village.
We want to find or design an algorithm to optimally swap land plots between us and other owners in order to best consolidate.

Is there such a tool in QGIS?


texts: QGIS displaying features in reversed order

texts:
I loaded in a .shp and went to inspect the attribute table, the features are listed as follows:

#  Attr_1  Attr_2
1       3       C
2       2       B
3       1       A  <--- hover shows "Feature ID: 0"

The table in ArcGIS however shows them in the order they are in their respective files:

#  Attr_1  Attr_2
1       1       A
2       2       B
3       3       C

I'm exploring QGIS so this was an interesting find. Is this normal?


texts: Seeking tool/algorithm for assigning code to Enumeration Areas (polygons) using serpentine numbering method?

texts:
One of the projects I am working on has asked me to assign a code (starting from 1 and counting up) to polygons that represent enumeration areas so the statistician can sample the areas. I have been asked that the assigning of the code uses a serpentine numbering method. The image below shows an example of what I am trying to achieve. You start in the North West part of the area and snake bank and forth assigning numbers.

Now I've searched various sites to see if anyone had developed a tool or even an algorithm that could be implemented that would achieve this. All I have found are numerous powerpoint slides and pdf files that simply say they have used serpentine numbering. I can only assume they have done this manually.

Does anyone know of an existing tool that would take a shapefile and number as described?

It seems to me that someone may have attempted this before but I have not been successful in finding a suitable reference as I don't usually get involved in census type data processing.


texts: How to add a new column in attribute table with attribute with indexes of another max-min ordered attribute field?

texts:
I've an attribute column named e.g. VALUE with different values I would like to make a column named e.g. ORDERsuch that the greatest value has the attribute 1 the second value 2 etc. etc. so order express the order from the biggest to the smallest value like this:

VALUE   ORDER
12000     1
8000      2
600       3
500       4
200       5

So far I've tried with QGIS 1.8 to order the VALUE column (that please note is not in the right order by default) from maxto minsimply clicking over the attribute name and than use $rownumin field calculator but this seems use the default rows numeration so I don't have the right order.

Please edit my question if you think you can explain it better!


texts: Where can I find cellmast data in the UK?

texts:
I am looking for a shapefile containing the locations of cell masts that provide phone coverage in the UK. On sites like Orange, O2 etc they have an interactive signal checker (https://explore.ee.co.uk/coverage-checker for example), but I haven't found a way to obtain the data in GIS form. Have rung several phone companies, but the customer services down there don't seem to understand what I am asking for!

Any help would be greatly appreciated, thanks.


texts: Issues Accessing ESRI File Geodatabase within Quantum GIS Desktop

texts:
I am trying to access ESRI file gdb from Quantum GIS. I am on Windows 7, 64-bit. I have found many useful threads on how to set up this functionality within Quantum GIS. The method I went for is described in this thread. I installed QGIS as part of the OSGeo4W package. I went for an 'Advanced' installation and added QGIS (v.1.8.0) + 'gdal-filegdb:OGR File GDB Driver' library. This has been installed with no issues. In QGIS, when I navigate to 'Add vector layer' window and choose File GDB as source type, I am able to preview contents of the file gdb in the following 'Select vector layers to add' window. Upon making my selection I can see the layer being displayed in table of contents, however  after a short while qgis is crashing. I have attached couple of screenshots of this process:

Does anybody have similar experience or else can suggest solution to this issue?
Do I need the ESRI File Geodatabase API installed as well to be able to access ESRI file GDB?

Many thanks,
Magda


texts:
Is it possible to reset or clear the map and plot new points? I have a search box that am using to query a geojson API.

After the first mapped result, the second search is plotted in addition which I do not want. What I want to achieve is, if I search, that result should be plotted not all subsequent searches I make.

Is there a way, to overwrite that first plotted result or remove it before plotting a new search result?

button.addEventListener('click', function(e) {
    inputVal = input.value;
    // console.log(inputVal);
    getPointLayer = new ol.layer.Vector({
        source: new ol.source.Vector({
            url: 'http://0.0.0.0:5000/api/v1/dumela/' + inputVal,
            format: new ol.format.GeoJSON()
        }),
        style: styleFunction
    });

    var layersOnMap = map.getLayers();
    // reset the map, clear all point
    if (layersOnMap) {
        map.removeLayer(getPointLayer);
    }

    console.log(layersOnMap.length);
    map.addLayer(getPointLayer);
});


texts: Update service based on enterprise database

texts:
So I have an enterprise database on SQL Server which contains a feature class (only points)

I add this feature class to a new ArcGIS Desktop project, so there is only this layer in the project.

I then go to File -> Publish as a service and choose to publish it on ArcGIS Server.

I then go to the service definition and check on "ArcGIS.com Map View" that the service was published correctly.

Then comes my problem :

I add a new point to the Feature Class in the database, I check that the new point was added to the Feature Class without issues and it displays correctly on ArcGIS Desktop when I add the layer to the project (taking it again from the enterprise database).

I then check again the service published on the server, but the modifications I did in the database AFTER publishing the service aren't displayed.

Why does the changes I made in the Feature Class stored in the database (which the service is based on) aren't synchronized on the published service ?

edit : the layer owner is the user "DBO" and the layer was taken from the database connected as "DBO" (through ArcMap UI). The publishing was made by a user made specifically for publishing services. The layer isn't versioned and the database is on arcgis server. When I published I didn't change the default settings.


texts: When I add DTED data files (*.dt1 and *.dt2) to my ArcGlobe application, they are displayed as solid gray boxes.  Do I need to do something special to the files and/or the ArcGlobe control to see these properly?


texts: ArcGIS JavaScript API 3.3 change the value of a field to lowercase in InfoTemplate

texts:
Is there a way to change the case of a text field in an InfoTemplate using ArcGIS JavaScript 3.3?

Let's say I am bringing in a feature layer, and it is a requirement that I ensure that the field value for {$gaugelid} is lowercase, as I am constructing a URL to preview that particular Gauge ID's Hydro-graph in the Info template.

The formatting for the Hydro-Graph URL is: https://water.weather.gov/resources/hydrographs/Lowercase(${gaugelid})_hg.png)

If it's not possible when constructing the InfoTemplate, is there a way to define that I want those values as lowercase in the definitionExpression or when I define the outFields?

var HydroURL = 'https://water.weather.gov/resources/hydrographs/${gaugelid}_hg.png';
var iframeHydro = "";
console.log(iframeHydro);
var content = "Location: ${location}
Status: ${status}
Observation:${observed}${units}
Observation Time: ${obstime}
"+iframeHydro;
var gaugetemplate = new InfoTemplate("${gaugelid}", content);

var currentGaugeData_lyr = new FeatureLayer(
    "https://idpgis.ncep.noaa.gov/arcgis/rest/services/NWS_Observations/ahps_riv_gauges/MapServer/0",
    {
      definitionExpression: "state = 'TX'",
      outFields: ["*"],
      infoTemplate: gaugetemplate
    }
  );
var gaugeSymbol = new SimpleMarkerSymbol();


texts: Python equivalent of spatstat's spatial point process functions

texts:
Is there a free Python package (i.e. PySAL?) that implements the standard summary statistics for spatial point processes, e.g. the F, G, J, and K functions?
I'm happy to use R (i.e. spatstat) if that's the best option, but I want to try out some other things I already have implemented in python. The implementations don't have to be nearly as exhaustive as R, just looking for something to get started.


texts:
I have a table osm_roads that has a column named geometry geometry geometry(LineString,900913).

If I want to filter a city to draw districts, I cannot do it by geographic coordinates?

What if in openJUMP or QGIS I want to select by bounding box, such a table as that?

Thank you in advance.


texts:
I have recently started using Geopandas and am having a problem with the performance of intersects. The intersect tests whether land parcels have a road within (or across the boundary).

Parcel data (parcels) is a polygon geodataframe and roads (roads) are as a line geodataframe. I attempt to find out if the parcel interacts with a road with the follwoing code:

road = roads.unary_union
parcels['road_intersection'] = parcels.intersects(road)

Here I generate get a column 'road_intersection' with True/False values.

with 985 parcel polygons and 20854 road lines, this takes over a minute to generate a result.

I have previously tried with

parcels['road_intersection'] = parcels.intersects(roads)

but this did not seem to test all roads with all parcels.

Would anyone be able to provide any tips for correct use or optimisation?


texts: Difference tool becomes unresponsive and causes QGIS (2.18.24) to crash

texts: I'm trying to find the difference between two layers (the first is a manitoba wide shapefile and the second is a layer that narrows down the water bodies of manitoba selected from another canadian waterbodies shapefile). Both layers have the same CRS (after modification). I'm trying to create a new layer that shows me only the land area of manitoba (minus the water bodies) but the difference tool keeps crashing on me. What exactly is going wrong?


texts: Leaflet dynamic legend for active layer controls?

texts:
I´m trying to achieve a map legend that shows the names of the active Layers and also adds new names to the list every time a new layer is checked.

Also it should only appear when at least one layer is active and disappear when no layer is checked.

Right now I´m pretty close to my goal I think but something is still not working correctly and I don´t know what it is.

var arr = new Array ();

var layerLegend = L.control({position: 'bottomleft'});

layerLegend.onAdd = function (map) {
var div = L.DomUtil.create('div', 'legend');
div.innerHTML +=
'Legende
'+test();
return div;
};

map.on('overlayadd', function (arr, eventLayer) {
        arr.push(eventLayer.name);
        layerLegend.addTo(map);
}
);

function test (arr) {
for (var i = 0; i < arr.length; i++) {
console.log(arr[i]);
}};

map.on('overlayremove', function (eventLayer) {
    map.removeControl(layerLegend);
}
);


texts: Where can I set the BufferColor property for labels in QGIS?

texts: The accepted answer to this question suggests that the buffer color can be set as a data defined label property. However, in the Layer Properties > Labels > Advanced tab, the only options I can see for the buffer are size and transparency. How can I set the BufferColor property?


texts: Impossible to launch QGIS after installation

texts: I deleted QGIS and I installed it again. Once I did that I receive the typical message saying "QGIS has been installed on your computer" but once I close the wizard I cannot find any icon to launch it, neither in the Windows start menu (I have windows 10), neither in the QGIS 3.10 folder. I've tried different alternatives (32bytes instead of 64, older versions, etc) but nothing works: it's being impossible to find a way to launch QGIS.


texts: How to draw radius of a circular point buffer using PostGIS

texts:
I have generated circular buffers of points in a shapefile. The generated circles and points share a common ID/key, say pt_id.

My question is

Is it possible to use PostGIS to generate a third layer of polylines (straight lines) representing the radii as follows?

I need this for better presentation.


texts: gdalnumeric.LoadFile() with temporary output from QGIS algorithm

texts:
I'm trying to access a previously created temporary raster output with gdalnumeric.LoadFile() to calculate some statistics. It seems to me that gdalnumeric.LoadFile() can only access a file through a file path and not with a QGIS temporary output.
I tried without success:
import gdalnumeric
path="C:\Users\...\xxx.gpkg"
raster = processing.run("gdal:rasterize", {'INPUT':path,'FIELD':'','BURN':1,'UNITS':1,'WIDTH':50,'HEIGHT':50,'EXTENT':'32351398.3019399,32593108.6666274,5366514.87416986,6045057.7298552 [EPSG:5652]','NODATA':0,'OPTIONS':'','DATA_TYPE':0,'INIT':None,'INVERT':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
raster_file = raster['OUTPUT']

Is there a way to access the temporary raster with gdalnumeric.LoadFile()? Or do I need to use a QGIS native algorithm to continue my analysis?


texts: Red x marks along lines after adding a length field?

texts:
I added a "length" field to a shapefile of lines called "roads".

QGIS puts red x's along the length of the road afterward.

Why is this, and can I make these go away?


texts: 'invalid parameters' on ArcGIS rest API used from Silverlight

texts:
I have an Silverlight application from which I'm trying to create new features to ArcGIS server. Every time I get error

{\"error\":{\"code\":400,\"message\":\"Unable to complete  operation.\",\"details\":[\"Invalid parameters\"]}}

My request code is like this:

System.Net.WebClient client = new System.Net.WebClient();
client.Headers["Content-type"] = "application/x-www-form-urlencoded";
client.Encoding = System.Text.Encoding.UTF8;
client.UploadStringCompleted += client_UploadStringCompleted;

client.UploadStringAsync(new Uri("http://.../FeatureServer/0/addFeatures?f=json&features"), "POST", JSONString);

If I copy the JSONString variable to a text editor, replace every escaped double quote mark (\") with double quote mark (") and use the ArcGIS servers /FeatureServer/0/addFeatures form, I am able to store the same features in that form to a featurelayer. That makes me thinking that the JSON data is well formed and valid for the featurelayer I'm using.

What could make the request response an error every time when called ArcGIS server from Silverlight application?

SOLVED:

I followed the instructions below and found a solution. First, the endpoint URL is like /.../FeatureServer/0/addFeatures, not like /.../FeatureServer/0/addFeatures?f=json?features

features=[{"geometry":{"x":0,"y":0}, "attributes":{"Id":"018-00","SV":"","Num":"4 A"}}]&f=json

Found the solution from here and by comparing the POST methods made from my application and the web form with Fiddler


texts: What could be use cases for adding an alpha band to a raster?

texts:
In QGIS I just used the tool "cliprasterbymasklayer", which produces a new raster layer as output and gives me the option to "create an alpha band". I read that the alpha band represents the transparency but I don't see any use case. Even if I don't choose this option I can make the new raster more transparent. If I check the option the new raster is lighter and has 0 as lowest value.

I was just wondering about a situation in which that could be practical. Any ideas or examples?


texts: ArcGIS Script Tool Input Parameters giving Error

texts:
I am using Desktop 10.5.1

I have a python script, works great when I copy/paste into the ArcMap terminal and replace the arcpy.GetParameterAsText function with a hard-coded string variable.

However when I put this script as a script tool inside a toolbox I get errors for the input values.

the code is below for the script tool.

import arcpy
arcpy.env.overwriteOutput = True

structures = arcpy.GetParameter[0]
Prim_Conductor = arcpy.GetParameter[1]
TableOutput = arcpy.GetParameter[2]
PrimConductor_FL = "Primary_Conductor FL"
structures_FL = "structures_lyr"
arcpy.MakeFeatureLayer_management(structures, structures_FL)

whereSQL = """ Damage_Cat = 'Major Damage' OR Damage_Cat = 'Minor Damage' """

arcpy.SelectLayerByAttribute_management(structures_FL, "NEW_SELECTION", whereSQL)

bufferOutput = r"\\path\to\shapefile.shp"
arcpy.Buffer_analysis(structures_FL, bufferOutput,"20 Feet", line_side="FULL", line_end_type="ROUND", dissolve_option="NONE", dissolve_field="", method="PLANAR")
arcpy.SelectLayerByAttribute_management(structures_FL, "CLEAR_SELECTION")
arcpy.MakeFeatureLayer_management(Prim_Conductor, PrimConductor_FL)
arcpy.SelectLayerByLocation_management(PrimConductor_FL, "INTERSECT", bufferOutput)

arcpy.TableToExcel_conversion(PrimConductor_FL, TableOutput, "NAME")

Parameter data types in the tool properties are:

[0]: Feature Class (it is a SDE featureclass)
[1]: Feature Class (SDE Feature Class)
[2]: File (output direction)- output for excel file.

The error I see is:

Traceback (most recent call last):
  File "Y:\Emergency Data\Electric V-Look Up\Script Tool\Script.py", line 7, in <module>
    structures = arcpy.GetParameter[0]
TypeError: 'function' object has no attribute '__getitem__'

Failed to execute (Script).


texts: Is there a method to tell which hemisphere a point belongs to, in WGS84 geoid, where the equator is an arbitrary ellipsis?

texts:
I use the WGS84 geoid (because I'm working with a GPS chip).

I define an ellipsis that splits the surface in two parts, lets name these 'A' and 'B'.

This ellipsis is defined by two points 'a' and 'b' , using their coordinates in WGS84. 'A' is the surface left of the vector 'ab', and therefore 'B' is right of vector 'ab'.

How do I tell if point 'c' belongs to 'A' or not? Any language will do, even pseudo code. I'll figure out how to translate to ATMega assembly later...

Many thanks.

PS: What I currently think of is switching origins so that my ellipsis is the new equator, then applying the transform on my point and finding out if it's North or South. But whereas I could find out how to do it in polar coordinates systems, I'm totally lost in geoids... maybe its not even possible?

PPS: Here is an example:

The points are:

50°21'37.9"N 3°05’01.8"E : A
50°21'36.3"N 3°05'04.2"E : B
50°21'40.2"N 3°05'05.9"E : C1, belongs to A
50°21'32.1"N 3°04'52.7"E : C2, does not belong to A


texts: Geoserver and Geofence - using custom rules

texts:
I've checked all the options of the rule creation in GeoFence Embedded (integrated with GeoServer) and GeoFence Standalone web application. The available options for building data access rules seem to be limited to: username, user's role, service name, request name, workspace and layer. The additional option is to limit the access with the geographical area.

Is there a way, using GeoServer and GeoFence, to make rules more specific, as in the following examples:

Restricting access to a group of features from a layer, which belong to the selected bigger layer (show only cities from a selected state)
Using CQL expressions in GeoFence rules with system variables (showing only features from a layer which have in one attribute the username of the logged-in, using GeoFence authenticator, user)


texts: PostGIS interpolation / triangulation options

texts:
I'm looking for an open source workflow to automate DEM construction.  We have a series of sites that fall within a LIDAR dataset.  We want to create site specific DEM's for each site, and we're looking to automate the process.

So far, we have automated:

data loading into PostGIS (LIDAR and site points)
creation of site boundaries (combination of st_buffer and st_envelope)

Now we're looking for options to interpolate the data subsets at each site and export them to Surfer7 grid files.

Currently, we're querying the data in the PostGIS database using QuantumGIS, exporting to csv files, then manually importing and gridding the data in Surfer7.  Hopefully we can automate this as well.

So for this use case, we'd like to triangulate our lidar data to a DEM.  Other potential cases we can think of involve different interpolation methods - so if there are options for inverse distance and kriging - we're interested!

This is very much a learning exercise for us - we're working in baby steps!


texts: Does anyone know how to display a netCDF in a 3D Scene Map in ArcGIS 10.2? something like the pretty 3d graphics at the UC I've got hundreds of thousands of points/cells with values for a stretch of time - and I'd like to display/play with symbolization from one file - though at worst I could display each year as from shapefile (cube symbol) and extrude one year atop another - though that sounds tedious.


texts: What are common functions GIS web app should have?

texts:
I'm a programmer making his first steps in the GIS world. In order to get to know this world I wrote some map based apps using WMS / OSM , OpenLayers and MS-SQL (non spatial).
Now I want to put it all together and do some more advanced coding (adding caching capabilities for example , thematic map...)

What are the most common functions in your opinion a GIS app should have ?

I get the obvious are measuring (line and area) probably something related to features info ... But beside that I have no idea what else.

I would like to use your experience for making a base app that will answer most basic GIS requirements (and will put together all the technologies I chose to practice with)


texts: Issues with displaying layer with the right extent and units

texts:
I am using QGIS 3.14 and I want to project a layer with a particular CRS, extent and units.
Starting situation: I have 2 shapefiles that are in different projections, have different extents and are displayed in different units.

The first layer (fishnet) has the (correct) projection of EPSG:32617 - UTM Zone 17N - Projected, extent in meters and units in meters.
The second layer (polygon) has a projection of EPSG:4326 - WGS 84 - Geographic, extent in decimal degrees and units in degrees.

When I both import them into QGIS, the layers are on top of each other.
Issue: Since the second shapefile is slightly tilted (due to different projection), I want to change its information to reflect that of the first shapefile. When I right-click on the layer in the Table of Content, click 'Set CRS' --> 'Set layer CRS' --> WGS 84 / UTM zone 17N, the shapefile disappears from field overview i.e. the shapefiles are no longer overlapping. When I right-click on the newly projected shapefile and click 'layer properties', I read that its CRS was indeed correctly changed as well as the units, but the extent is still decimal degrees. When I right-click the same layer and click 'zoom to layer', it appears somewhere completely different.
Any suggestions on how I could solve this issue?


texts: Writing test suite for QGIS standalone application?

texts:
I have a rather simple QGIS standalone Python app (2.8.2, Windows 7/64), with some published methods to create different types of geometric primitives (rectangles, circles).

class MyApp(QMainWindow, UI_MyApp)
    def __init__(self, ...)
    def circle(self, x, y, radius)
        ...
        layer.triggerRepaint()
    ...

def run_tests()
    self.circle(...)
    self.circle(...)
    ...

def main()
    ...
    app = QApplication(sys.argv)
    map = MyApp()
    map.show()
    map.raise_()
    app.exec_()
    ...

Now I want to run some automated tests against these methods, in this example circle().

It doesn't matter if I include the test calls inside the class itself (as shown above, e.g. connected with a button in the UI) or kept in a separate module: in both cases all graphics are not drawn until run_tests() finished.

How can I implement a (separate) test module which can call methods of MyApp() application with immediate response after each call and thus acting like an API?


texts: Import multiple WorldImages in a single GeoServer store

texts:
I downloaded tiles from Copernicus EU DEM.  Each tile is made of four files :

$ ls *E30N30*
eu_dem_v11_E30N30.TFw
eu_dem_v11_E30N30.TIF.aux.xml
eu_dem_v11_E30N30.TIF
eu_dem_v11_E30N30.TIF.ovr

I guess these are WorldImage files ?  Now I am trying to create a store from them in GeoServer.  However the WorldImage store type only seems to accept single WorldImages, and the ImageMosaic store type fails to load data with

Could not list layers for this store, an error occurred retrieving them:
Failed to create reader from file:///data/source/Copernicus_EU_DEM_v1.1 and hints
Hints: REPOSITORY = org.geoserver.catalog.CatalogRepository@433abb55
_EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@13eb5bc5[Running, pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0]
System defaults:
FILTER_FACTORY = FilterFactoryImpl
FEATURE_FACTORY = org.geotools.feature.LenientFeatureFactoryImpl@2081b1a2
LENIENT_DATUM_SHIFT = true
FORCE_AXIS_ORDER_HONORING = http
COMPARISON_TOLERANCE = 1.0E-8
GRID_COVERAGE_FACTORY = GridCoverageFactory
TILE_ENCODING = null
STYLE_FACTORY = StyleFactoryImpl
FORCE_LONGITUDE_FIRST_AXIS_ORDER = true

(I had a similar issue with another dataset earlier because the dataset directory was missing write access for tomcat user.  This is not the same case here, the permissions are right.)

How am I supposed to load these tiles ?


texts:
I noticed some inconsistency on geographies between ST_Project and ST_Segmentize in postgis when dealing with trying to segmentize over very long distances (around 18000-19000 km)

Below are two plpgsql functions that should give about the same result but they don't. I tried it using different versions of Postgis (up to 2.1.8).

For both functions I use a start point lat and long, an azimuth, a total length and a segment length.

segmentize_line_manually: I project the starting point at incremented distances (segment length) with the given azimuth using ST_Project.
segmentize_line_auto: I project the starting point at the total distance with the given azimuth, then I call ST_Segmentize with the segment length.

Around 18000-19000 km I notice a rather important difference between the two lines. In my opinion, the segmentize_line_manually gives the correct result.

CREATE OR REPLACE FUNCTION segmentize_line_manually(startLat DECIMAL,

startLong DECIMAL, azimuth DECIMAL, max_distance_meter INT, interval_meter INT)
RETURNS geometry AS
$$
DECLARE
    line_points geometry[];
    startPoint geography := ST_MakePoint(startLong, startLat);
    distance_meter INT := 0;
    i INT := 1;
BEGIN
    WHILE distance_meter < max_distance_meter LOOP
        distance_meter := distance_meter + interval_meter;
        line_points[i] = geometry(ST_Project(startPoint, distance_meter, azimuth));
        i := i + 1;
    END LOOP;

    RETURN ST_MakeLine(line_points);
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;

CREATE OR REPLACE FUNCTION segmentize_line_auto(startLat DECIMAL, startLong DECIMAL, azimuth DECIMAL, max_distance_meter INT, interval_meter INT)
RETURNS geometry AS
$$
DECLARE
    line_points geometry[];
    azimuth_line geometry;
BEGIN
    line_points[1] := geography(ST_MakePoint(startLong, startLat));
    line_points[2] := ST_Project(line_points[1], max_distance_meter, azimuth);

    azimuth_line := ST_MakeLine(line_points);
    RETURN st_segmentize(geography(azimuth_line), interval_meter);
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;

I call both methods using the following parameters:

SELECT segmentize_line_****(0, 0, 0.628973508651991, 19000000, 50000);

The results below. In red, the line created using ST_Project, in blue with ST_Segmentize:

Is there a problem with my code, or is it really an issue with ST_Segmentize?


texts: Creating raster layer from numpy array using pyqgis?

texts:
I am working on a plugin for Qgis to calculate spatial Kernel density maps. I have all the calculations working, all I am missing is a way to turn a Numpy Array, with density values into a multiband raster layer.

Do I have to create a geotiff on a temp file using Gdal and then load it?

Or is there a direct way to create the layer from data in memory?

if so, how to do it?


texts: Seeking workaround for completely_within spatial join in ArcGIS Desktop?

texts:
I tried to perform a 1-to-many spatial join between polygons (target feature) and lines (join features) to join the lines that are completely covered by a polygon to that polygon. (the polygon is actually a result of buffering the lines with a dissolve option)
I used the option 'completely contains', but this returned a background processing error. After sending parts of my data to Esri Inc, I got the confirmation that there is indeed a bug they need to fix (#NIM086228).

Anyway, my project needs to be completed. In search of a work-around, I already tried to turn the lines into polygones (small buffer around it), but it returns the same error. When trying the contains_clementini (or actually all contains and within options) suddenly give me an error even before running it. The error is 'Relationship invalid for selected layers'.

Here is what I need to do:

  I have (pipe)lines that I want to 'merge' when they have common
  characteristics. So I create a buffer and let it dissolve when the
  value of those characteristics is equal. A polygon layer is created,
  but there is not link back to the underlying lines that are the basis
  of the polygon. In order to have these ID's, I wanted to perform the
  spatial join with match option 'completely contains

Has someone any experience with the last type of invalid error messages or another workaround for me?


texts: should ST_Centroid and ST_Transform be applied in a particular order?

texts:
I am a complete newbie using GIS, and particularly using postgresql-postgis.

I have a table with administrative boundaries where limits are multipolygons defined in utm (EPSG:23030).

I want to calculate both centroids and also minimum-bounding-circle centre and radius, for each row of the table.

To my surprise, the coordinates I got are different if I transform them before/after calculating the centroid:

postgis=# SELECT gid, name, ST_AsText(ST_Transform(ST_Centroid(the_geom),4326)) as centr1_trans2, ST_AsText(ST_Centroid(ST_Transform(the_geom,4326))) as trans1_centr2 from adminboundaries_23030 where name='Madrid';

 gid  |  name  |              centr1_trans2               |              trans1_centr2
------+--------+------------------------------------------+-----------------------------------------
 5642 | Madrid | POINT(-3.6962798422859 40.4655298819556) | POINT(-3.6963536465966 40.465587716731)
(1 row)

If I do it the other way round (conversion from a table with latlon 4326 to utm 23030), these are my results. 6-8 m of difference:

postgis=# SELECT gid, name, ST_AsText(ST_Transform(ST_Centroid(the_geom),23030)) as centr1_trans2, ST_AsText(ST_Centroid(ST_Transform(the_geom,23030))) as trans1_centr2 from admboundaries_4326 where name='Madrid';

 gid  |  name  |              centr1_trans2               |              trans1_centr2
------+--------+------------------------------------------+------------------------------------------
 5642 | Madrid | POINT(441074.289262773 4479875.27865048) | POINT(441080.495696979 4479868.80970194)
(1 row)

I know this is not a huge difference, but is there any easy explanation for it?

Which one would be the most correct way of calculate the centroids before displaying them in a web map?


texts: Turning addresses into GPS coordinates into delivery route?

texts:
I have a list of locations (prospects I need to meet) I have to drive to, and I'd like to use either Google Maps or Open Street View to build an optimized route.

Are there tools that can take a CSV list of locations and generate a driving route?

Thank you.


texts:
I'm very new to raster manipulation, and I'm getting a nodata value of -3.4028234663852886e+38 when importing USGS NED 1/3 arc-second n41w074 1 x 1 degree IMG using rasterio on Python.

What is it with this number? I understand that it's some minimum bound of a 32bit float, but why is this number so common in raster data? Why not just, say 0, or -999 or something less extreme?
Is it a problem? How to deal with it? Google provides little answers about why. There is one thread that kind of deals with this but it focuses on how to solve it using QGIS (not my case) or preprocessing with GDAL. I've seen a few posts of using a mask, which I understand inherits from numpy masks and conceals the numbers for the purposes of calculations, etc.

I'm trying to merge two rasters and I'm getting the wrong results (as is clear in the bottom figure) so I'm thinking it might be related to this.

What should I do?


texts: Tool in Toolbar is not visible in Customise Commans tab

texts:
I have added a tool in a toolbar in python addin wizard. And while opening arcmap I am getting the error :---

c:\arcgis\ArcToolbox\Toolboxes\Space Time Pattern Mining Tools.pyt
^

SyntaxError: invalid syntax

and tool is not visible in command tab of customize window. Please help me how to solve this problem.


texts:
Possible Duplicate:
  GIS QUERY on point and buffer

All,

I have a PostGIS/Postgresql DB that has one table in it. One with point geometries and the other that represents buffer as polygons. i would like to updted the flag field in the same table when ever mine point location geometry is changed. if the point is outside the buffer polygon then its value should be 0, or flag value should be 1 when it is inside the buffer polygon.

mine table is as follows

gid flag    point(geom)         polygon(geom)

1   2   the point Geometry  the polygon Geometry
2   1   the point Geometry  the polygon Geometry
3   2   the point Geometry  the polygon Geometry
4   1   the point Geometry  the polygon Geometry

manish


texts: Equivalent to ArcPy's PointGeometry in ArcObjects?

texts:
I'm trying to select a county feature from a nationwide counties layer using Geoprocessor's managed SelectLayerByLocation coclass:

internal static void GetCountyName(IPoint centerpoint, out string countyName)
{
    Geoprocessor gp = new Geoprocessor();
    SelectLayerByLocation selectByLocation = new SelectLayerByLocation();
    selectByLocation.in_layer = countiesLayer; //IFeatureLayer in mxd
    selectByLocation.select_features = centerpoint; //cause of error
    try
        {
            gp.Execute(selectByLocation, null);
            //countyName will be retrieved from value of selected feature
        }
    catch (Exception)
        {
            string message = "";
            for (int i = 0; i < gp.MessageCount; i++)
                message += gp.GetMessage(i) + "\n";
            MessageBox.Show(message);
        }
}

The MessageBox showed me this:

I found here that the "select_features" parameter requires either a FeatureLayer or FeatureClass:

But I know that ArcPy has PointGeometry specifically so that you don't have to create a stored feature class just to use it as a geoprocessing parameter:

  In many geoprocessing workflows, you may need to run a specific operation using coordinate and geometry information but don't necessarily want to go through the process of creating a new (temporary) feature class, populating the feature class with cursors, using the feature class, then deleting the temporary feature class. Geometry objects can be used instead for both input and output to make geoprocessing easier.

Is there an option like this for ArcObjects?


texts:
After about 10 hours of dissolving, QGIS gives me this message :-(

Any ideas, what this means and maybe how to correct?


texts: How to remove abnormalities from the Green vegetation fraction

texts:
Study area was clipped by using the shape file from Band 3 and Band 4. First, Band 3 and Band 4 were converted to at-satellite brightness temperature by using the built-in function in semi-automated tool in QGIS (semi-automatic classification plug-in).
NDVI was calculated by b4-b3/b4+b3 in landsat 5 images, the range of NDVI was 0.0968299 - 0.433846.
Green Vegetation fraction were also needed to examine the role of vegetation in mix land cover that was derived from the equation:

  GVF = NDVI - NDVImin / NDVImax - NDVImin

but the problem here was as you would see the range of the fraction is -0.631684 - -0.0011024) that is due to some abnormalities in the data that don't know how to remove it. Is there something wrong in "NDVI" values as there is no data values (no data values = -3.4020 + 38) by seeing the properties of image in QGIS.
I have taken the NDVImin and NDVImax as 0.0968299 and 0.433846 respectively but literature has shown that the values 0.05 and 0.70 as NDVI threshold values that were derived by global Advanced Very High Resolution Radiometer (AVHRR) 10-day composite data (Zeng et al. 2000).


texts: QGIS python processing.runalg Clip returns None

texts:
I'm writing a script in pyQGIS, which would emulate measurement of a line by "walking" along it with a compass with constant length. I have a problem with the Clip algorythm from processing module, namely it returns None, regardless of whether I add layers as inputs or selected features. What do I do wrong? I believe, that should the clip result in an empty layer, it would still be there, just empty? Below is the piece of code with problems,

while d.measureLine(QgsPoint(coorFin[0],coorFin[1]),QgsPoint(coor[0],coor[1]))>step: # If the distance between current and last point is more than chosen step
    # Create buffer polygon around a point
    bufCirc = point.geometry().buffer(step,10) # do a buffer from existing point
    buforek = QgsFeature()
    buforek.setGeometry(bufCirc)

# Create layer and put the polygon inside
    bufferLayer = QgsVectorLayer("Polygon", "temporary_polygons", "memory")
    pz = bufferLayer.dataProvider()
    bufferLayer.setCrs(my_crs)
    pz.addFeatures([buforek])
    bufferLayer.commitChanges()

# select the buffer
    buforSel = bufferLayer.setSelectedFeatures([buforek])

# Clip coastline with the buffer
    liniaWarstwa = processing.runandload("qgis:clip", selectedCoastLine, buforSel, "memory:liniaWarstwa")
    print liniaWarstwa
    layer = QgsMapLayerRegistry.instance().mapLayersByName("memory:liniaWarstwa")[0]

The print prints "none", QgsMapLayerRegistry says "index out of range".


texts:
CartoDB LayerActions Toggle

Is there a way to make one button that shows a sublayer on one click, then hides on the next click. So one button per sublayer which turns the layer on and turns the layer off. Similar to the functionality built into the CartoDB UI.

For example how could I reprogram a toggle button by using some other piece of code, or could I expand the function in the HTML file by using a flag value (0 for an ON state which could be show(), and another value such as 1 for an OFF state which would be hide())

This is my code below:

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/themes   /css/cartodb.css" />

Novant Health
All


texts: No data found in bounds using rioxarray to clip GeoTiff to shapefile


texts: AttributeError: 'int' object has no attribute 'setName' (or insertChildNode)

texts:
I am trying to make group and fill it with layers in pyQGIS.

Piece of my code:

from PyQt4 import QtCore, QtGui
from PyQt4.QtCore import QFileInfo

from qgis.core import *
from qgis.utils import iface

group = iface.legendInterface().addGroup( 'abc')
group.setName("Group X")

But I got only error AttributeError: 'int' object has no attribute 'setName'

May I import some other class or where is problem?


texts:
I replaced some gs-wps (geoScripting files from C:\Program Files (x86)\Boundless\OpenGeo\jetty\webapps\geoserver\WEB-INF\lib) and stopped the geoserver and again started it. After that it giving me the error.
HTTP ERROR: 503
Problem accessing /geoserver/web/. Reason:

Service Unavailable


texts: Enabling Input Features for Clipping Geometry on a Raster in ModelBuilder using parameters?

texts: I have a model to summarize tree heights and density using LiDAR.  It utilizes model parameters for inputs.  I can't enable the Clip option to Input Features for Clipping Geometry because the output extent is derived from a parameter.  How do I assign a default value to the model to enable the clip option while still having the model run from parameter inputs?


texts: Exporting well data from IHS Kingdom in a format that plays nice with ArcGIS

texts:
A coworker of mine has a bunch of well data (including georeferencing and attribute tables) in IHS Kingdom.  We need to convert this data into a format that works with ArcGIS (e.g. .shp files).  Unfortunately, it looks like the only options for exporting are as follows:

.wtv2 - Tobin Well

.asc - Landmark Well 2.x/3.x
.asc - Landmark Well OWX

.wel - Charisma

.gds - GDSII

.dat - Well Production Data

As far as I can tell, none of these are compatible, and each of them are extremely poorly documented online.  As is Kingdom IHS, for that matter.

The only solution we've managed so far would be to export the map as a large PDF before georectifying and manually digitizing each of the wells in ArcGIS.  However, there are too many wells to make this a viable option.


texts: Concatenate multiple string / text fields some of which may have NULL values ArcGIS

texts:
I am trying to repeat the behavior of this code from the ESRI website.

I want to replicate the output example in ALL_Type field, so that it concatenates a field but it deals (ignores) any that have NULL values.

Unfortunately this example no longer works at 10.2.2 because all example of NULL return a "None" text. Instead I just want a blank to show.

And here is the code.

   "*args" allows this routine to accept any number of field values.
# the values are passed as a Python tuple, essentially a
# non-editable list
#
def concat(*args):

  # Initialize the return value to an empty string,
  # then set the separator character
  #
  retval = ""
  sep = "_"

  # For each value passed in...
  #
  for t in args:
    # Convert to a string (this catches any numbers),
    # then remove leading and trailing blanks
    #
    s = str(t).strip()

    # Add the field value to the return value, using the separator
    # defined above
    #
    if s <> '':
      retval += sep + s

  # Strip of any leading separators before returning the value
  #
  return retval.lstrip(sep)


texts: How Do I Add Custom "Region" labels with GeoJSON polygon style datapoints, similar to country names

texts:
I have a Simple CRS map I've set up in Leaflet, and I'd like to label portions of the map with just text labels, similarly to how on a typical World Map you'd see all the continents, then zooming in it shows countries, then states, essentially having text labels only for certain zoom levels.

Is there a way to create these labels and automatically center them with GeoJSON coordinates for polygons around the regions I would like to label?


texts:
I'm trying to get familiar with Folium library in Python. I wanted to make a web-map of Turkey. Framework uses EPSG:3857 Web Mercator CRS as default. I tried to use approximate coordinates at the centroid of Turkey in EPSG:4326 CRS which is lat:39.0 and lon:35.0 by specifying EPSG:4326 as CRS while giving the parameters but output map pointed out somewhere out of projection area near Pacific Ocean. I looked at epsg.io to confirm that I'm not confusing CRS but could not figure what is the problem myself. Below are my two examples:

1- crs; EPSG:3857 coordinates;(lat:39, lon:35)

Centers at Syria, unsurprising to have misalignment since CRS is Web-Mercator and coordinates at Geographic WGS84

webMap = folium.Map(location=(39,35),
                    zoom_start=5,
                    crs='EPSG3857') #which is default

# visualize
webMap

2- crs; EPSG:4326 coordinates;(lat:39, lon:35)

Centers at out of projection area

webMap = folium.Map(location=(39,35),
                    zoom_start=10,
                    crs='EPSG4326')

 visualize

webMap

Can anyone point out what I am missing?


texts: st_warp a star_proxy object with multiple files

texts:
I have a s stars_proxy object with 12 attributes in files, i want to reproject each files in a new crs and resolution stars_proxy object, i use this command in R :
r <- stars::st_warp(src = s, dest = d, use_gdal = TRUE, method = "bilinear")
But in r i have just 1 attribute, not the 12 attributes file reprojected ?


texts: Merging multiple MrSid files in ArcGIS for Desktop?

texts:
I've got multiple sid files that make up my county of interest.

I need to get all the sid files into one and STILL have the sid file type after merging.

Anyone know how to do this in ArcMap 10.2.2?


texts: Combining two layer symbologies into single legend item using QGIS?

texts:
In QGIS 2.18.13 on Windows 10, I have symbolized a polygon with two layers on top of either other; the bottom layer is a solid fill color, and the top layer is a BDiagonal (hash).  The two layers together produce a colored polygon with a hashed overlay in the map.

In the print composer, is there a way to combine the two legend entries that result from these two layers into a single element, such that the symbol for that polygon is shown as a single colored box with the hashing over top of it?

Right now I just have the two separate layers with their separate symbology.  I need them combined.

Alternatively, maybe there is a way to produce this same symbology in a single layer on the map prior to the composer, so that this isn't an issue.  I just haven't found a way to do this.


texts: Approach to smoothly transitioning a large amount of GeoJSON features' style

texts:
I'm working on an interactive Leaflet map with about 1000 GeoJSON features. When I click on a one of the menu items on right, I use a function to style the feature according the values of the active property:

           let toChangeTo = getColor(layer.feature.properties[propName], activeKey, subdivisionDefaultStyle);
           layer.setStyle(toChangeTo);

This is the CSS I use so that the features transition:

  .leaflet-interactive{
    transition: fill .75s;
   }

I've made sure my GeoJSON files are as small as possible. I'm wondering if there's a way to change the styles of the features successively, rather than all at once, so there isn't so much lag. Can anyone recommend an approach to make this animation run at a faster FPS?

Below is a demonstration:


texts: How to do a WFS version 2.0.0 GetFeature request with OL3?

texts:
How to obtain information about the feature WFS layer in different versions of WFS 1.0.0 / 1.1.0 / 2.0.0 ?

The OL documentation says

  Feature format for reading and writing data in the WFS format. By default, supports WFS version 1.1.0. You can pass a GML format as option if you want to read a WFS that contains GML2 (WFS 1.0.0). Also see ol.format.GMLBase which is used by this format.

Ok, example WFS - GetFeature is performed in the version 1.1.0

var f = ol.format.ogc.filter;
var wfs = new ol.format.WFS();
wfs.writeGetFeature({
    srsName: 'EPSG:3857',
    featureNS: 'http://openstreemap.org',
    featurePrefix: 'osm',
    featureTypes: ['water_areas'],
    outputFormat: 'application/json',
    filter: f.and(
        f.like('name', 'Mississippi*'),
        f.equalTo('waterway', 'riverbank')
    )
});

But it always forms a request in WFS version 1.1.0.
How to form GetFeature request in WFS versions 2.0.0 using writeGetFeature function?


texts: Setting up PostGIS database and opening it in QGIS on Windows?

texts:
I have tried to follow the manuals for setting up a PostGIS database to be accessed by users with QGIS, but at the moment I cannot understand enough to progress?

My situation; I am a trainee responsible for creating a PostGIS database for the digitising of green objects (such as trees) and running queries - like linking a street list, with image files, to the attribute table.

I have "shapefile skeletons" (currently empty with no objects, just the attribute frame" where the current and future georeferenced objects (a different shapefile for each type) will be pasted/created). I trust I can just import those skeleton shapefiles into the database?

The situation;

We work in a separate building with our own LAN and internet connection with a backup server for every file. So we have local file access and files that we share over an internet connection

The setup of the database has been done by an administrator via remote access, using the latest .exe execution files of PostGIS and PostgreSQL (postgresql-9.6.2-3-windows-x64 & postgis-bundle-pg96x64-setup-2.3.2-1) for my 64-bit machine. The administrator provided me with a text file with login information. The training manuals here are pretty confusing since they deal with Ubuntu.

I am stuck as I do not know where to use the login information that I have written down? I assume I need to use pgAdmin?


texts: How to enable Map Tip display text in qgis and qgiscloud?

texts: I'm using the free version of qgiscloud and qgis 2.4. I noticed that if I'll use the HTML option in "Map Tip display text", the map in qgiscloud will only show the title of the layer when I hover/click but it's working if I used the Field option. I would like to use the HTML option since I want to display/style multiple values based from an expression from different attributes of my layer, instead of using the Field option which will only display the selected field.


texts: Loading raster on PostgreSQL/PostGIS Linux server not working

texts:
I used the following raster2pgsql.exe to upload an SRTM mosaic to my Linux database server (working fine for vector data):

raster2pgsql.exe -I -t auto e:\data\Sinai\SRTM90V41b\SRTM41M90.tif SRTM90c | psql -d Sinai-U postgis -h 179.19.1.27 -p 5432

The process is carried out smoothly, however, when I open the DB Manager in QGIS, I'm unable to preview the Public.SRTM90c raster, and if I try to add it to QGIS, I get this message:

Invalid Layer: Raster layer Cannot instantiate the 'postgresraster' data provider

Any clues?


texts: Changing pixel value in raster layer in QGIS using Python?

texts:
I'd like to change the value of a pixel in a QGIS raster layer.  But the following Python code cannot do it.

r_layer = iface.activeLayer()
provider = r_layer.dataProvider()
from osgeo import gdal
raster = gdal.Open(str(provider.dataSourceUri()))
np_array = raster.ReadAsArray()

np_array[1,2] = 5

raster.GetRasterBand(1).WriteArray(np_array, 0, 0)
raster.GetRasterBand(1).FlushCache()


texts: ArcGIS Online won't display weather widgets

texts:
Trying to add a custom pop-up using real time weather underground personal weather station data feeds like:
https://www.wunderground.com/personal-weather-station/dashboard?ID=KFLBUSHN3

When I add these in a data field or as a url link via an html file it works fine in ArcGIS Desktop. They want a shared ArcGIS Online app instead. Doesn't seem to work in there. i can get the url in a pop-up, but no display of the data feed happens..


texts: I am trying to have markers in a KML to follow the direction of the view in Google Earth

texts: I am trying to represent the flow direction of streams in a KML file to be used in Google Earth. As complex line symbols are not supported by the KML format, I have opted for a solution combining lines and point symbols. In ArcMap 10.1, I have populated a new field in my stream layer with the direction of the lines in degrees. I have used Feature to point to generate points. Since KML don't support point rotation either, I have made eight categories based on the line direction in degrees (North, Northeast, East, etc.) and I associated corresponding arrow symbols. So far so good, but when I opened the KML file in Google Earth, I realized that the symbols are right just as long as the view points North up. If the view changes direction (e.g. South up), my symbols do not adapt. Here (https://productforums.google.com/forum/#!topic/earth/I668BTY7AP4), someone suggests a heading to the icon style (of 56 degrees in the example). What I would like is to have my symbols follow Google Earths view , or in other words, North, rather than an orientation in degrees.


texts: Arguments for opening a GIS area in an important agricultural-related bank?

texts:
the managers of the bank where I work at have been considering creating a "GIS-unit" that would help with economic-related research / business intelligence. They asked me to raise a list of reasons of why/how such group would benefit the bank businesses. So, my friends, I pass this question to you all. Which arguments can you think of for making the case for the bank to create such group? Just a note, this bank is strong in agricultural-related investments.

TL/DR: How can GIS benefit an agricultural-bank businesses???


texts: Set my Home location on Google Maps, if it does not have an official/registered address?

texts:
The answer to "Adding new address to Google Maps?" does not include a case for non-business addresses. I want to add my home address, which is not a register address. Ideally, I should have been able to add a home (or other private-use) location, from a point on the map. This functionality does not seem to exist. I do not want an official location, visible to everybody that would be suitable for a business. I want to set a private location and name its address as I wish (because there is an address, used locally).

Can I do any of this? Any ideas?

Update-rephrase because "unclear":

In Google Maps, under My Places there are some "Labeled" personal locations, one can set, which will be known on all Google products. Among them there is "Home" and "Work" and also the choice to add others. Before there was only the choice to set the address of those "Labeled" places, in some text-boxes (plus some other staff for the place, always in text). Now, there is also a location pointer on the map, which you can move to specific coordinates. The issue has been resolved, just after this post (the next day) and there is an answer already; that is how "unclear" the question was.

[Observe that this update does not do anything to "clarify your specific problem or add additional details to highlight exactly what you need," which is stated, clearly enough, in the original question: (1) "I should have been able to add a home (or other private-use) location, from a point on the map."; (2) "I want to set a private location and name its address as I wish (because there is an address, used locally)". The only imprecision, in my opinion, is that I do not stress that the locally used address is a real address, despite not being registered and, since my "Home" is a private location, it might have been possible to set the address itself --and maybe search for it! (This second part has not been resolved or answered. I am not sure I like giving such ideas to Google, but someone would eventually...)]