texts:
I'm looking for a feature in QGIS but I am not sure if it exists.
I don't know the name of such feature, but I call it multiple selection sets. It exists in TransCAD and I find it very useful. It basically allows one to save the selected features under a name, and have multiple sets that one can compose (union, intersection, difference, etc).
I am under the impression that QGIS allows only to save the selected set as a new file, which is not what I am looking for.
Does it exists in QGIS? How can one use it? If not, is there a work-around?


texts: Creating iterative Select by Location, Summary Statistics tool using ModelBuilder?

texts:
I do a lot of proximity analysis using "select by location" and would like to cut down on the effort spent each time doing this!

I have put together a model using ModelBuilder that:

Iterates over features in a feature layer (source layer);
Selects features in a separate target layer based on the set parameters; and
Uses summary statistics to summarize over a certain field.

Where I'm getting caught up is in making this as much like a standard tool as possible.

I have two issues I'm running into:

I've set the "statistics field" of the Summary Statistics tool as a model parameter. However, when I go to run this tool through ArcToolbox, I set the target layer, but the statistics field isn't populated with the list of potential fields from the target layer. The drop down menu remains blank.
Is there a way to make all the target layer fields populate similarly to how they would if you were using the Summary Statistics tool alone and set an input layer?
Get Field Value gets the value of the first row of the specified field, but I do not know how to tell it to get the value of the field of what I just summarized on.

So, in Summary Statistics, if I told it to get the SUM of a field called POPTOTAL, then I want the Get Field Value tool to grab the value from the new field called "SUM_POPTOT" in the output statistics.
There should only be one row of data, so I don't need to worry about getting the correct row. I have "Field" set as a model parameter currently, but I don't think it needs to be.

I've posted this thread in the ArcGIS forum, but I haven't received any feedback. I'm still new to ModelBuilder...I have some background in programming, so Python script would be fine (or maybe even better).


texts: I'm looking for a tidal range GIS layer. Anything global that shows either mean tidal range, maximum tidal range, or maximum tidal height. I have been trawling the web for tide related GIS layers and they don't seem to be available. If any one knows of a source, please let me know!


texts: Moving labels manually using expression in QGIS?

texts: I have many point shapefiles that are too close together.  I need to label 4 lines of data for each point. I have written the expression that shows the 4 lines of data and now the labels are overposting on top of eachother. How can I manually move the labels to another location. I have used the data defined coordinate X,Y options to do this and told the program to use the Expression, please see the attached image, but it does not activate the move label button in the label toolbar. The label toolbar only activates if I choose 1 attribute at a time from a list.


texts: Prevent stretching across the entire map when connecting coordinates into path from Asia to North America

texts:
I am connecting coordinates around the world in a specific order. When moving from Asia to North America the line makes a connection that stretches across the entire map instead of bleeding over the edge of the map and reconnecting in North America. How can I avoid it?
I have three layers:

a basemap (ESRI Gray)
CSV data
output of the "Points to Path" geoalgorithm

I found this answer, but it seems complex and even if that is the solution to my issue as well I have been unsuccessful at implementing it in my case: QGIS Vector layer stretched across entire map


texts: How to increase FUSION-LTK lidar processing productivity?

texts:
I am processing Lidar data with FUSION-LTK.  My study area contains 5 blocks of 1000 las tiles per block.  Each block is broken into approximately 4 subblocks of 250 las tiles.  I am calculating grid metrics and strata metrics to assess vegetation structure.  I wrote simple batch files to process entire blocks (e.g. block 1 subblock1; block 1 subblock2;...) using the command line tools in FUSION-LTK.

The problem I am encountering is that I am only able to process one subblock at a time.  When I try to simultaneously start processing on another block, the current session of FUSION-LTK stops the processing and moves to the new block.  How can I improve the processing workflow so that more than one subblock are processed simultaneously, or at least enhance the multiprocessing capabilities of LTK?


texts: Clip raster using a shapefile of 5500 polygons

texts: I have a nightlights raster file of India, and i want to clip it using the block shapefile for India. There are 5500 blocks. But I can't find a unique attribute(except for serial no.) to clip the raster with, I tried doing it using ModelBuilder in ArcMap 10, but if I put the serial no. in "group by fields" there is no output(the rasters are generated but it is a blank image when i put it in arcmap). If I put some non-unique atribute, the model runs for sometime then shows error.


texts: Implementing spatial tables in a database

texts:
Is there a OGC specification on how to implement spatial data storage in a database?

I was looking at the tinydb for spatial data storage, but I was wondering what is a standard way to store spatial data.


texts: OpenLayers 3 bug when rendering GeoServer WFS layers?

texts:
I'm having issues rendering WFS layers with OpenLayers 3. If I follow the example here and zoom out zoom out zoom out, my pc will soon freeze as the amount of data is huge. Don't even want to try on a phone.

If I set MaxFeatures to e.g. 500, the example works well. However I then alter the code to run with MY GeoServer, and I receive:

Assertion failed. See https://openlayers.org/en/v4.6.4/doc/errors/#36 for details.

The text for error #36 is:

Unknown SRS type. Expected "name" or "EPSG"

My GeoServer is returning JSON as requested. I check the JSON sent from my GeoServer:

... "type": "FeatureCollection", "totalFeatures": 0,
    "features": [{"type": "Feature", "id": "plats.66", "geometry": {
        "type": "Polygon",…},
    ],
    "crs":{"type":"EPSG","properties":{"code":"3857"}}}

So, my GeoServer returns "type": "EPSG", however the OL WFS request parser seems to only accept "type":"name", this in spite of the error message.

My GeoServer is v2.5.3, the WFS request is v1.1.0 as per the OL example.

Is this a bug, if so is there a workaround?


texts:
I'm interested in generating a dbf from a feature class in a gdb. I am also interested in selecting only 6 out of the many fields, which is why I added steps to keep only wanted fields. I am using python on ArcMap and ArcCatalog 10.3.1. Here is the code:

fc = r"C:\Users\migrate\Desktop\SFBG_GIN_Test\gdb\SFBG_GIN_Test.gdb\PlantCenter_Test"
dbfLocation = r"C:\Users\migrate\Desktop\SFBG_GIN_Test\dbf"
dbfOutputName = "Plants_GIS"
dbfPath = dbfLocation + '\\' + dbfOutputName + '.dbf'

if arcpy.Exists(dbfPath): # This line exists to overwrite any existing dbfs
    dm.Delete(dbfPath)

arcpy.TableToTable_conversion (fc, dbfLocation, dbfOutputName)

# Get all fields in dbf and remove unwanted fields
fields = arcpy.ListFields("fc_temp")
keepFields = ['GIS_2016_08_25', 'SourceAcce', 'AccessionC', 'SectionNam', 'Longitude', 'Latitude']
dropFields = [x.name for x in fields if x.name not in keepFields]
dm.DeleteField (dbfPath, dropFields)
print "DBF created"

I get the following error:


texts: How to do a Kernel Density with "population field" in r?

texts:
In ArcGIS, I can generate a kernel density map and weigh it by population field. For example, I can use city's population as weighting factor, so that more populated city has higher influence to the surrounding area.

I am trying to do the same in r. Right now I am using the package Splancs. However, I can only generate a kernel density map by treating every point the same. In other words, I cannot give more weights to certain points (e.g., a larger city).

Please let me know how to do it in r, and which package should I use.


texts:
I am trying to get an user input for getting a parcel name and pan the map to that selection. What am I doing wrong in the code below?

import arcpy

#Msg box to get parcel Name
parcel_Name = str(arcpy.GetParameterAsText(0))
arcpy.AddMessage ("Loading parcel " + parcel_Name)

try:

    mxd = arcpy.mapping.MapDocument('CURRENT')
    df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
    ssparcel = arcpy.mapping.ListLayers(mapDoc, "ssparcel", df)[0]

    whereClause = "Name ='%s'" % parcel_Name
    arcpy.SelectLayerByAttribute_management(ssparcel, "NEW_SELECTION", whereClause)
    df.extent = df.panToExtent(ssparcel.getSelectedExtent())
        #df.scale *= 1.5
    arcpy.RefreshActiveView()

except:

    arcpy.AddWarning ("Error")

...

Removed the "try" and got the following error:

Executing: Script "0002"
Start Time: Thu Jun 05 11:32:14 2014
Running script Script...
IndentationError: unexpected indent (abcd.py, line 9)
Failed to execute (Script).
Failed at Thu Jun 05 11:32:14 2014 (Elapsed Time: 0.00 seconds)


texts:
Is there any reporting tool that can be used to build custom PDF reports with maps from OpenLayers?
I tried iReports, but i was not able to create such report. Any suggestions?


texts: Overlay briefly has incorrect position in OpenLayers

texts:
I have an OpenLayers map with a number of points on the map. Some of those points are labeled via an ol.Overlay. Every 10 seconds, I call to the server and ask if the location of that point has changed, and the labels are re-created.

In the first image below, the application is paused right after the locations have been refreshed and the new label is created, but is not shown yet (the old label will be shut off just prior to showing the new label). The inspector is highlighting the hidden but placed label.

Point #170 hasn't changed position, but initially, the label created for it is situated 13px to the right of where it should be.

I hide the old label, and show the new label (still appearing in the incorrect place) :

But, after a fraction of a second of being shown, the new label moves to its correct place:

This creates an annoying flickering appearance in some browsers (Chrome is fine, but in the others it is more apparent).

Is this a result of OpenLayers, or is it more likely to be an issue in my code?

The problem seems to resolve itself after this line (83093) in ol-debug.js:

this.dispatchEvent(
  new ol.MapEvent(ol.MapEventType.POSTRENDER, this, frameState));


texts: Unable to load Shapefile using pgShapeloader / shp2pgsql-gui (pgAdmin3's PostGIS Shapefile Import/Export Tool)

texts:
I am trying to load a shapefile (http://dds.cr.usgs.gov/pub/data/nationalatlas/statesp020_nt00032.tar.gz) into PostGIS 2.0.1 using pgAdmin3's PostGIS Shapefile Import/Export Tool.

However its not allowing me to import the shapefile, showing the error messages below. Any idea what happen?

Output from PostGIS Shapefile Import/Export Tool

==============================
Importing with configuration: statesp020, staging, geom, C:\gisdata\statesp020_nt00032\statesp020, mode=c, dump=1, simple=0, geography=0, index=0, shape=1, srid=4269
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Failed SQL begins: "SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;
BEGIN;
CREATE TABLE "staging"."statesp020" (gid serial,
"area" float8,
"perimeter" float8,
"statesp020" float8,
"state" varchar(20),
"state_fips" varchar(2),
"order_adm" int4,
"month_adm""
Failed in pgui_exec(): ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 15: SELECT AddGeometryColumn('staging','statesp020','geom','4269...
                ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Shapefile import failed.


texts: How do I patch libgeos without losing QGIS?

texts:
How can I create a custom libgeos to replace the one in my Ubuntu system?

Basically, I'm trying to expose some of GEOS functions (orientedDistance()) to  other packages as explained in Directed Hausdorff distance in PostGIS?.

I guess I can trying compiling GEOS from source. But the problem is that, even if I succeed in figuring out a patch for exposing the function, I don't know how to replace the current libgeos in my system without breaking things.

To test, I tried removing the Ubuntu provided package with

sudo apt-get remove libgeos-3.5.1 libgeos-c1v5 -s

Then from the result below, I can see that most of the GIS packages I need will be gone:

...
0 upgraded, 0 newly installed, 35 to remove and 3 not upgraded.
Remv qgis-plugin-grass [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
Remv qgis-provider-grass [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
Remv qgis [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
Remv python-qgis [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
Remv python-qgis-common [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
Remv gdal-bin [2.2.1+dfsg-2~xenial0]
Remv grass-core [7.2.1-1~xenial4]
Remv postgresql-9.6-pgrouting [2.4.1-1~exp1.pgdg16.04+1]
Remv postgresql-9.6-postgis-2.3 [2.3.3+dfsg-1.pgdg16.04+1]
Remv python3-gdal [2.2.1+dfsg-2~xenial0]
Remv postgis [2.3.3+dfsg-1.pgdg16.04+1]
Remv python3-geopandas [0.1.1-3]
Remv python3-fiona [1.7.6-1~xenial1]
Remv libgdal1i [1.11.3+dfsg-3build2]
Remv qgis-providers [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis]
...
Remv libgdal20 [2.2.1+dfsg-2~xenial0] [libqgis-core2.99.0:amd64 libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv libqgispython2.99.0 [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis] [libqgis-core2.99.0:amd64 libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv libqgis-core2.99.0 [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis] [libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv python3-shapely [1.5.13-1build1] [libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv libgeos-c1v5 [3.5.1-3~xenial0] [libsqlite3-mod-spatialite:amd64 libspatialite7:amd64 libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv libgeos-3.5.1 [3.5.1-3~xenial0] [libsqlite3-mod-spatialite:amd64 libspatialite7:amd64 libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv liblas-c3 [1.8.1-4~xenial1] [libsqlite3-mod-spatialite:amd64 libspatialite7:amd64 libqgis-customwidgets:amd64 liblas3:amd64 ]
Remv liblas3 [1.8.1-4~xenial1] [libsqlite3-mod-spatialite:amd64 libspatialite7:amd64 libqgis-customwidgets:amd64 ]
Remv libqgis-customwidgets [1:2.99.0+git20170820+f867b65+24xenial-ubuntugis] [libsqlite3-mod-spatialite:amd64 libspatialite7:amd64 ]
Remv libspatialite7 [4.3.0a-5] [libsqlite3-mod-spatialite:amd64 ]
Remv libsqlite3-mod-spatialite [4.3.0a-5]

I plan to also compile PostGIS from source, so that's not a problem.

But how can I keep other packages such as qgis, gdal etc. intact? Is this possible?


texts: Perhaps I'm missing something obvious - but is there no built-in spell check tool in QGIS? I know spell checking attribute data may be a waste of time; but I would expect a tool for layouts, just to ensure labels etc have a basic sort of check would be available, am I missing something?


texts: Uncaught TypeError: Cannot read property 'split' of null

texts:
I am getting this error when i display "name" on the map from geojson for "point" entities.

var pointStyle = (function () {
  var style1 = new ol.style.Style({
    image: new ol.style.Icon({
      scale: 1,
      src: imgPath + 'location.png'
    }),
    text: new ol.style.Text({
      text: null,
      scale: 1,
      fill: new ol.style.Fill({
        color: '#000000'
      }),
      stroke: new ol.style.Stroke({
        color: '#FFFF99',
        width: 3.5
      })
    })
  });
  var styles = [style1];
  return function (feature, resolution) {
    style1.getText().setText(feature.get("name"));
    return styles;
  };
})()

When i run my angular 2 application, map is rendered with point entities (some point entities with "name"). All the entities are not displayed with "name". And this is when i get this error.

When i zoom in , then all the point entities "name" is loaded. Then i don't get this error.

How to solve this ?


texts: How to create tests for QGIS' processing models

texts:
Short version: Is there a way to create automated tests inside my own processing models?

Extended version:

I am just starting with QGIS' processing to create some complex model which makes quite a lot of changes and I would like to know if there is a way to test my own model while launching it, so if the tests fail, I can get a notification and I will know that something went wrong before I realize too late and I have to lose any changes I made after executing the model with unexpected results.

Example:

I am doing several joins and spatial joins. I learnt the hard way that after triggering several times my model, every time I was joining layers, new fields were added with a _1 appended to them in order not to overwrite the previously created fields by mistake when I repeated the model.
It would be great to have some test that either checked all the layer's fields or counted its number and, if they weren't as expected, triggered an alarm/dialog box or created a log file.

For the little I know, I assume the way to achieve this is by creating a Python script and add it to my model as one step more, but I don't have any clue on how should that script be in order to interact with the model.


texts: Programmatically changing raster symbology after reclassification using ArcPy

texts:
I have created a code in Python that I will use for a custom tool in ArcGIS, the code will create a Eucledian Distance of 6000m around a point and Reclassify the distance from 0-2000, 2000-4000, 4000-6000. The code works but what I was wondering was, is there a bit of code I can add on to make the 0-2000 area red, the 2000-4000 area to be orange and the 4000-6000 area to be yellow. As the code stands the colours are random i guess (blue, green etc.). I have been looking through a few books but could not find a solution.

from arcpy import env
from arcpy.sa import *

# Set environment settings
arcpy.env.workspace = "Z:\\GIS TEST\\Select_by_Location"
arcpy.env.overwriteOutput=1

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Local variables:
Breakout_Location = "Breakout_Location.shp"
maxDistance = "6000"
Output_direction_raster = "Z:\\GIS TEST\\Codes\\Output"

# Set Extents
arcpy.env.extent = arcpy.Extent(31652.290929,28102.389486,269652.513962,212966.447186)

if arcpy.Exists(Output_direction_raster):
    arcpy.Delete_management(Output_direction_raster)
Euc_Dist1=EucDistance(Breakout_Location, 6000, 400, Output_direction_raster)

# Process: Reclassify
Reclass_Output_direction_raster = Reclassify(Euc_Dist1, "Value", RemapRange([[0,2000,1],[2000,4000,2],[4000,6000,3]]), "DATA")

arcpy.SetParameterAsText(0, Breakout_Location)

Reclass_Output_direction_raster.save("Z:\\GIS TEST\\Codes\\Output")


texts: Mouse pick point snap to isnappingfeedback location

texts:
I have created one BaseTool for placing point on polygon feature. mousemove event showing mouse move time vertex & edges & endPoints but picked point snapped at vertex.

 public override void OnMouseMove(int Button, int Shift, int X, int Y)
        {
            IMxDocument pMxdoc = m_application.Document as IMxDocument;
            IMap pMap = pMxdoc.FocusMap;
            pPosition = pMxdoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y);
            if (pPosition != null)
            {

                ISnappingResult snapResult = null;
                //Try to snap the current position
                snapResult = pSnapper.Snap(pPosition);
                pSnappingFeedback.Update(snapResult, 0);
                ISnapAgent agent = pSnappingEnv as ISnapAgent;
                if (snapResult != null)
                {
                    pPosition = snapResult.Location;
                }
                else
                {
                    IPoint point = pMxdoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y);
                    var activeView = pMxdoc.ActiveView;
                    point = pPosition;
                    var Chamged = psnapEnv.SnapPoint(point);
                }
            }
            // TODO:  Add Tool1.OnMouseMove implementation
        }

 public override void OnMouseDown(int Button, int Shift, int X, int Y)
        {
            IMxDocument pMxdoc = m_application.Document as IMxDocument;
            IMap pMap = pMxdoc.FocusMap;
            ESRI.ArcGIS.Carto.IActiveView pview = pMxdoc.FocusMap as IActiveView;
            pMap.ClearSelection();
            pPosition = pMxdoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y);
            if (pPosition != null)
            {
                for (int i = 0; i < pMap.LayerCount; i++)
                {
                    ILayer pLayer = pMap.Layer[i];
                    IFeatureLayer pFlayer = pLayer as IFeatureLayer;
                    IFeatureClass pFclass = pFlayer.FeatureClass;
                    if (pLayer.Name == "Edit Comment")
                    {
                        ESRI.ArcGIS.Geometry.Point pntclass = pMxdoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y)as ESRI.ArcGIS.Geometry.Point;
                        CreateFeature(pFclass, pPosition, pFlayer);
                    }
                }
            }
            // TODO:  Add Tool1.OnMouseDown implementation
        }


texts: In the QGIS error message "Feature (99) has invalid geometry" does the 99 refer to the OBJECTID?

texts:
QGIS 3.4.2-Madeira running simplify on a polygon z shapefile with 1075 rows. "Feature (99) has invalid geometry" not sure how to find this feature. If I edit and delete the row with OBJECTID = 99 and rerun simplify the same error presents but lists Feature(98) as invalid.
Interestingly running Geometry Tools\ Check Validity returns no errors and no invalid output, but if I run the fix geometries tool and then run simplify on the Fixed geometries layer then simplify works fine.
Key question I have is how to identify (and fix) the problem geometry for "Feature 99". Does the 99 refer to the OBJECTID?


texts: QGIS raster size is different from size in file?

texts:
I created a GeoTIFF file with size 3600x3600, I verified that the file I have on harddrive is indeed 3600x3600. But the raster size that QGIS displays is 4554x2277.
This makes the raster sum wrong, (407Million instead of the correct 509Million). Why does QGIS show a different raster size here? Why the particular number 4554?
Script I used for creating the raster file.
from osgeo import gdal
from osgeo import osr
import numpy
import math

fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)

xdim = 3600
ydim = 3600

dst_filename = "test.tif"
dst_ds = driver.Create(dst_filename, xsize=xdim, ysize=ydim,bands=1, eType=gdal.GDT_Float32)

x0 = -180
y0 = -90
x_per_pixel = 360.0 / xdim
y_per_pixel = 180.0 / ydim

dst_ds.SetGeoTransform([x0, x_per_pixel, 0, y0, 0, y_per_pixel])

srs = osr.SpatialReference()
#srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("WGS84")
proj_str = srs.ExportToWkt()
print('proj: ',proj_str)
dst_ds.SetProjection(proj_str)

raster = numpy.zeros((xdim, ydim), dtype=numpy.float)
earth_a =  6378137 / 1000
earth_b = 6356752.3142 / 1000
for i in range(xdim):
    for j in range(ydim):
        lat = y0 + y_per_pixel * (j + 0.5)
        a = 2 * math.pi * earth_a / xdim
        b = math.pi * earth_b / ydim
        raster[j][i] = math.cos(math.radians(lat)) * a * b

dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
dst_ds = None

What QGIS shows in Properties->Information of the raster layer
Name    test
Path    I:\gis\worlddiv\py\test.tif
CRS EPSG:4326 - WGS 84 - Geographic
Extent  -180.0000000000000000,-90.0126558050849894 : 180.0253116101699788,90.0000000000000000
Unit    degrees
Width   4554
Height  2277
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description
Compression
Band 1
STATISTICS_APPROXIMATE=YES
STATISTICS_MAXIMUM=61.75239944458
STATISTICS_MEAN=38.018121292831
STATISTICS_MINIMUM=0.026944568380713
STATISTICS_STDDEV=19.920540978552
STATISTICS_VALID_PERCENT=100
More information
Dimensions  X: 4554 Y: 2277 Bands: 1
Origin  -180,90
Pixel Size  0.07905694150420948529,-0.07905694150420948529


texts: How to get the NDVI value and the coordinates of each pixel within a geometry?

texts:
I would like to know how to extract the NDVI values from a list of points and export them in a file in CSV format with the following structure: LAT / LONG / NDVI.
I managed to get to a certain point in the code and I couldn't move forward.
var collection = sent2
.filterDate('2015-06-30', '2016-06-30')
.filterBounds(geometry)
.sort('CLOUDY_PIXEL_PERCENTAGE')
.select('B2','B3','B4','B8');

// Image list
var imagesList = collection.toList(collection.size());

// Best image
var image = ee.Image(imagesList.get(6));

// NDVI
var red = image.select('B4');
var nir = image.select('B8');

var ndvi = image.expression('(nir - red) / (nir + red)', {
  nir: nir,
  red: red
}).select(['B8'],['ndvi']);

// Getting the projection of images
var proj = image.projection();

// Lat Long image
var latLong = ee.Image.pixelLonLat().reproject(proj);

// Getting a list of lat and long coordinates
var coords = latLong.select(['longitude', 'latitude'])
                      .reduceRegion({
                        reducer: ee.Reducer.toList(),
                        geometry: geometry,
                        scale: 10
                      });

// Putting lat and long information in separate lists
var lat = ee.List(coords.get('latitude'));
var long = ee.List(coords.get('longitude'));

// Creating a list of coordinates and converting them into points
var points = long.zip(lat).map(function (latLongList) {
  return ee.Geometry.Point(latLongList);
});


texts: Changing display properties of a layer using ArcGIS Server REST API Export Map function

texts:
I have a quick question I hope someone can help with.

Using the ArcGIS Server (10) REST API's Export Map function can I specify the boundary type, colour, fill colour and other properties in the features that are in the image?

I know I can easily do this if I use client-side rendering in Flex or Silverlight, and I have previously done this server-side in WMS with SLD, but so far I haven't seen this in a map image retrieved from the REST API.

There is some documentation on the Export Map function but I couldn't find what I was looking for.

Any info would be much appreciated.


texts: How to add row number based on distinct attribute value in QGIS (version 2.18.13)?

texts:
I have several layers (shapefiles) that have data in groups. Each group has a distinct group name. I am trying to number each item in the layer based on its distinct group name.

So, for each group I want to add a serial number for each object. For example group "broadst" = 1....7, group "crew" = 1......6 and so on.

Manually this can be done by sorting and selecting each group and adding a sequential number using rownumber.  But I have a large dataset and it will take a long time. Is there a plugin or would it be possible to use the field calculator using some form distinct filter based on group name?


texts: Is there a QGIS plugin that reclassify raster symbology based on current extent?

texts:
I am working with an elevation raster. I have assigned it color ramp symbology based on elevation. It looks like as the image shown below currently when zoomed to the layer.

But when I zoom in to start working on the map, the symbology is still the same; hence the color difference is not stark enough to differentiate between low points and high points.

My question is that; is there a plugin which can, with the help of one button press, reset my symbology based on my current extent so that the highest point in that area and the lowest point have entirely different colors.

In ArcMap, one can do that easily with the reset legend elevation option for TINs and Terrains as shown in the second picture below.

I am looking for something similar in QGIS for rasters. If there is no plugin already, can this kind of plugin be created using Python ?


texts: Why are the new vertex editing protocols not working properly in QGIS 3.16?

texts: I am using QGIS 3.16 to create new shapefiles from existing raster maps (p 225 in the 3.16 User Guide). My new polygons work OK when incorporated into other projects, but I cannot edit them to correct errors in their vertices. I am trying to use the protocol on page 485 in the 3.16 User Guide to (1) toggle editing, (2) left click an existing incorrect vertex once, (3) move the cursor to an improved position and (4) left click a second time, hoping to move the faulty vertex. Most of the time, I am able to select a vertex, but thereafter nothing happens at all, and the misplaced vertex remains in its original wrong position. On a few very rare occasions I have seen changes, perhaps after some unintended random change in my procedure.


texts: Solid is invalid : PolyhedralSurface (shell) 0 is invalid: self intersects


texts: How to import laser survey data to build raster. QGIS3.10


texts: Calculating distances with lat/long in decimals and WGS 84

texts:
I need to calculate the distance between points of two different datasets, once with locations of towns, the other with locations of dams. Both are located in different countries in Africa.

The more I read the more confused I keep getting. I'm using the Python ecosystem (Geopandas, Shapely, Fiona), but my question is basic enough that a general answer is helpful (if you provide some code it would be an added benefit!)

A) The first dataset is a .shp file with locations of towns as Points. This one is nice and provides the crs as epsg=4326 (which I understand is the code for WGS 84):

import geopandas as gpd
towns = gpd.read_file('towns.shp')) #Uses fiona to load
print(towns.crs)
print(towns.geometry.head(3))
#{'init': 'epsg:4326'}
#0    POINT (8.877318824300001 9.93427297769)
#1    POINT (9.163896418389999 9.47532028526)

B) The second one has locations of dams. It came in an excel file with latitude and longitude with decimal places. It didn't give any information on the crs.

#Omitting the excel loading of `dams`
print(dams.geometry.head(2))
print(dams.crs) # empty
#487    POINT (8.97333333333 9.76472222222)
#488    POINT (4.55305555556 8.44277777778)

This is my question

What's the right way to measure the distance between these two points? (Let's assume for now they are both in WGS 84)
I can calculate the distance between all towns and the first element of the dams:

# Distance function uses Shapely
print(towns.geometry.distance(dams.geometry.iloc[1]))
#0       0.194849
#1       0.346508
#2       1.046174

Is this just calculating the Euclidian distance and hence very inaccurate?
What do y'all do for this workflow? Should I transform the crs of the points to something that works for all of Africa and then take the Geopandas/Shapely distance function? Or would it be easier to keep the lat/lon (or WGS 84) and use a Haversine formula (or similar)? To me this would break a bit the benefit of using Geopandas.


texts:
I am reading in a shapefile of Metro areas in the US and then select the metro area of interest.
# Read shape file

px_shp = 'data/map/az/cb_2018_us_cbsa_500k.shp'

geodata = gpd.read_file(px_shp)
geodata.dtypes

# Select Phoenix-Mesa-Scottsdale, AZ

df_geodata = geodata[geodata['NAME'].str.contains('Phoenix-Mesa-Scottsdale, AZ')]
df_geodata.plot(figsize=(10, 6))
plt.title("Phoenix-Mesa-Scottsdale, AZ")

Now, I'd like to define the cities separately by passing the coordinates of the polygon. This means I'd have the find the boundaries of Phoenix, Mesa and Scottsdale cities. However, I cannot find a way/tool/API to obtain the polygon boundaries/geometries and pass them in.
I have looked at http://nominatim.openstreetmap.org/ which returns a collection of multipolygon in each of the cities.
[{"place_id":233116526,"licence":"Data © OpenStreetMap contributors, ODbL 1.0. https://osm.org/copyright","osm_type":"relation","osm_id":110815,"boundingbox":["33.277624","33.513376","-111.8939822","-111.5806478"],"lat":"33.4151117","lon":"-111.8314773","display_name":"Mesa, Maricopa County, Arizona, United States of America","class":"boundary","type":"administrative","importance":0.644252319667828,"icon":"https://nominatim.openstreetmap.org/images/mapicons/poi_boundary_administrative.p.20.png","geojson":{"type":"MultiPolygon","coordinates":[[[[-111.598158,33.45082],[-111.598195,33.45094],[-111.59823,33.451023],[-111.598257,33.451071],[-111.5983,33.45111],[-111.598342,33.451135],[-111.601778,33.451257],[-111.602976,33.451299],[-111.60446,33.451399],[-111.605305,33.451383],[-111.606666,33.4513881],

However, I am looking for aggregate of polygons within the multipolygon list - (POLYGON ((6065941.393835935 2104148.464510527.....)))
Ideally, I'd like a tool or an API that returns the polygon geometry based on name of the place.


texts: How can I calculate the shortest path between 2 points and insert the line result as a feature in a layer?

texts:
I try to make a script to find the shortest path between two bridges in a network. The goal is to get after each calculation a line that I insert as an entity in a predefined layer.
For now I do well to make calculations but when I change the coordinates to make another calculation nothing appears even if I heard the new entity that is inserted in my table.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

#building the graph**
vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

#coordinates of the Start point and endPoint**
pStart = QgsPoint(4.820590,45.838684)
pStop = QgsPoint(4.822459,45.837933)
tiedPoints = director.makeGraph(builder, [pStart, pStop])
graph = builder.graph()

#Calculation of the shortest path**
tStart = tiedPoints[0]
tStop = tiedPoints[1]

idStart = graph.findVertex(tStart)
idStop = graph.findVertex(tStop)

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)

if tree[idStop] == -1:
  print "Path not found"
else:
  p = []
  curPos = idStop
  while curPos != idStart:
    pnt = graph.vertex(graph.arc(tree[curPos]).inVertex()).point()

    p.append(pnt)

    curPos = graph.arc(tree[curPos]).outVertex()

  p.append(tStart)

  # Object generated after the calculation a qgsRubberBand**
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.red)

 #I build here a polyline based on the points in the p[] array
geom = QgsGeometry.fromPolyline(p)

for pnt in p:
    rb.addPoint(pnt)

 # I add here my result in a memory layer called cable and i want to do it for every feature i get after transforming the qgsrubberband object

v_layer = None
for layer in QgsMapLayerRegistry.instance().mapLayers().values():
   if layer.name() == 'cable':
        v_layer = layer
if v_layer is None:
    v_layer = QgsVectorLayer("LineString", "cable", "memory")
    v_layer.addAttribute(QgsField("id", QVariant.String))
    v_layer.addAttribute(QgsField("type", QVariant.String))
    QgsMapLayerRegistry.instance().addMapLayers([v_layer])

pr = v_layer.dataProvider()
v_layer.startEditing()

seg = QgsFeature()
seg.setGeometry(geom)
pr.addFeatures( [seg] )

v_layer.commitChanges()


texts:
I am creating a project with QGIS 2.18 to be exported to QField.
The Survey form is as follow (3 columns):

Customer name        - checked on site?  - Comment
Water meter number   - checked on site?  - Comment
Type of water meter  - checked on site?  - Comment
Tariff code          - checked on site?  - Comment
Diameter             - checked on site?  - Comment
etc. etc.

As one can sees, there will be several fields with the same properties ("checked on site" as form value relation, "comment" as a text).

Is it possible to duplicate some fields structure keeping the properties (widget, length, type, etc.) ?


texts: Creating ellipse with azimuth according to projected North not to geodetic North with ArcGIS Tool Table to Ellipse?

texts:
Table to Ellipse allows the user to create ellipses according to self-defined parameters. One of those parameters is azimuth. As written in the tool help the azimuth refers to the angle 0°-360° clockwise starting with 0° in North. I derived the parameters (coordinate of centre, widths and azimuth) from a projected map (UTM_34N) and used the corresponding projection information (*.prj) as the optional parameter Spatial Reference. As output I got a slightly rotated Ellipse see picture below.

After playing around with the tool I realised that the slight rotation corresponds to the difference between geodetic and projected North.

Do you see a possibility to create (not draw) an ellipse with the azimuth according to the projected North not to the geodetic North?

ArcGIS solutions are preferred but QGIS or other are very welcome too.

In the picture below you can see an ellipse, created as described above with the azimuth = 0°.


texts: I have put a buffer around points from a point data set. But some buffers overlap each other. Now I would like to remove  all the buffers (from the data set) which overlap another buffer. Does somebody have a suggestion, I would like to use the modelbuilder as I dont know how to use python.


texts: ModelBuilder model to classify Koppen Climate classification

texts:
I am trying  to see which area in the United States falls into which category of Koppen's climate zones.
Right now, I am working on B climates, using Raster Calculation.
The attached is a screenshot of my ModelBuilder.
PRISM_tmean is a raster data of average annual temperature in US, and PRISM_ppt is a total of annual average precipitation.
What I am struggling with is that now I want to use Con tool to say "if 70% of annual precipitation falls in summer half months (from April to September), use the Summervalue.
How can I use Con tool to do this?
I have raster data of average precipitation for all months as well.


texts: overpass query: select area/polygon with certain node in it

texts:
I want to select all polygons with

  tourism=camp_site and tourism=caravan_site without the key
  'sanitary_dump_station'

but where a node is located inside the boundary with the tag

  amenity=sanitary_dump_station

I tried different approaches and ended up with this:

 [out:json][bbox:{{bbox}}][timeout:120];

node ["amenity"="sanitary_dump_station"] ->.a;
.a is_in ->.b;
way(pivot.b)->.x;
way.x["tourism"~"caravan_site|camp_site"][!"sanitary_dump_station"]->.z;
(.z; .a; );
out geom;
>;
out skel qt;

In the results I find some campings without the key 'sanitary_dump_station', eg.: https://www.openstreetmap.org/way/116540519#map=19/50.22837/1.61213, while there is a node within its boundaries (https://overpass-turbo.eu/s/GRl)


texts: OpenLayers 2.12 + click handler + get the name of the vector

texts:
Basically,

I add several vector layers (KML content got from jQuery GET) in this way:

var layer1Vector= new OpenLayers.Layer.Vector('activeInfo');
layer1Vector.addFeatures(featuresLayer1Vector);
mp.addLayer(layer1Vector);

It works quite well.

Instead of adding the

selectControlKml = new OpenLayers.Control.SelectFeature(
    [layer1, layer2, ***], {
        clickout: true, toggle: false,
        multiple: false, hover: false,
        toggleKey: 'ctrlKey',
        multipleKey: 'shiftKey'
    }
);

I knew that it's possible to add a "general" click handler event

OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {
    defaultHandlerOptions:{
        'single': true,
        'double': false,
        'pixelTolerance': 5,
        'stopSingle': true,
        'stopDouble': false,
        'delay': 0
    },
    initialize: function(options) {
        this.handlerOptions = OpenLayers.Util.extend({}, this.defaultHandlerOptions);
        OpenLayers.Control.prototype.initialize.apply(this, arguments);
        this.handler = new OpenLayers.Handler.Click(
            this, {'click': this.trigger},
            this.handlerOptions
        );
    },
    trigger: function(e) {

        if (e.target._featureId) {
            //var feature = vectorLayer.getFeatureById(e.target._featureId);
            alert('click on fea' + e.target._featureId)
            //$(e.target).css('fill', '#000000');
            //$(e.target).css('cursor', 'wait');
            //$("div#info").append("You just clicked on " + feature.id + "
");
        }
***

If I click on a vector layer I have to open different custom popups (qTip2). The code:

if (e.target._featureId) {

works, but how can I get the layer name based on the e object? I analyzed all the methods/properties, but without having success.

E.g.
One of my layers was defined as:

var layer1Vector= new OpenLayers.Layer.Vector('activeInfo');

I need to get 'activeInfo'.

Thanks a lot!


texts: Creating line connecting polygons/points based on their location using QGIS

texts:
I have polygon layer which looks like this (here is small cutout):

What I need is to connect all polygons with line in somehow "logical" order - it means to go from first to last polygon (eventually point as I can create centroids) and create line connecting two nearest polygons BUT with one condition - there can't be closed polygon (it means I can't draw lines between three points so they create triangle) and there can be several branches - to make it clear below is the image how it should looks like.
What I get so far:

in QGIS I can't use Points to Path as I have no order field (all the polygons has same ID = 1)
in ArcMap I am able to use two function Generate Near Table and then XY to Line which connects at least some polygons together -->

finally I need to manually draw lines for the rest -->

Is there any automatic way to solve this using QGIS?
Maybe adding order field to layer with values based on location or something like that. I would ideally create an automatic model in the end as I deal with lots of these layers.


texts: Openlayers error: Uncaught TypeError: t.getLayerStatesArray is not a function

texts:
I am getting the following error in the yellow background below. Can anyone explain why i get this?
I am not using any bundler i include the ol.js  javascript file in the footer

  Group.js:209 Uncaught TypeError: t.getLayerStatesArray is not a function
      at Group.js:209
      at e.forEach (Collection.js:137)
      at e.getLayerStatesArray (Group.js:208)
      at e.renderFrame_ (PluggableMap.js:1198)
      at e. (PluggableMap.js:191)

        <script>
    var wmsSource = new ol.source.ImageWMS({
          url: 'http://webip/geoserver/wms',
          params: {'LAYERS': 'sweb:Streets'},
          serverType: 'geoserver',
          crossOrigin: 'anonymous'
          });

          var wmsLayer = new ol.source.Image({
            source: wmsSource
                });

                var baseLayer = new ol.layer.Tile({
            source: new ol.source.OSM()
                });

      var map = new ol.Map({
        target: 'map',
        layers: [baseLayer,wmsLayer],
        view: new ol.View({
          center: ol.proj.fromLonLat([33.21441650390625, 35.03449433167977]),
          zoom: 10
        })
      });

    map.on('click', function(evt){
        console.log(ol.proj.transform(evt.coordinate, 'EPSG:3857', 'EPSG:4326'));
    });
    </script>


texts: Saving edit widget setting in style or layer definition file?

texts:
I am using QGIS 2.4 on Windows 7 installed by the OSGeo4W-installer.

I want to save the "edit widget"-settings in a style to use them after loading in a different vector layer.

The problem: If I change the "edit widget"-setting in the layer properties (category "Fields") and save them in a style (QGIS layer style) or a layer definition file, the settings are not maintained when I load the style at a different vector layer (shp) or I open the layer definition file once again. It seems that the "edit widget"-settings are only be saved in the qgis-project.


texts: Sum habitat type within region using MapInfo?

texts:
I have two layers; one shows admin regions and the other shows habitats within those regions with columns for specific attributes and the corresponding area in hectares.

I need to work out the total area for each habitat attribute and I'd like to split this by region from the first layer.

Is there a simple way to do this in MapInfo?


texts: Weld individual line segments into one LineString using Shapely

texts:
I'm using Shapely in python and I'm given a MultiLineString with a bunch of Linestring objects. I can guarantee that all the LineString objects are simple lines with only 2 vertices and that they are all part of one single line (no branches).

I want to "connect-the-dots" and create a single LineString. Do I need to write a recursive welding method for this or is there a faster way?


texts: How to make gdalbuildvrt available in processing modeller?

texts:
I would like to make a model in QGIS that builds a virtualraster from several raster files. Unfortunately the gdal function is not available on the modeller. So, my question is: how to make this function available? It is necessary to construct a python script? If so how can I indicate the several raster files?
Thanks.


texts: Spatialite: how to add column with reference key?

texts:
I'm quite new to Spatialite but I need to build a small geographical database. I created a table and now I need to add a new column to it and to associate to this column a foreign key. Is it possible in any way?

I found out the code for SQLite, but not for Spatialite:

alter table Table1 add NewColumn integer references ExternalTable (id) on delete cascade

I've tried something that sounds good to me:

ALTER TABLE Routes
ADD COLUMN category INTEGER,
CONSTRAINT fk_routes_cat
    FOREIGN KEY (category)
    REFERENCES Vocabulary (id);

But Spatialite-gui gave an error:

  SQL error: "near "FOREIGN": syntax error"


texts: Speeding up my Python code for reading from Excel, performing geospatial analysis and writing to Excel files

texts:
I have written a Python snippet that reads lat and long stored in an Excel file. Converts them to a point which is then used to perform multiple geospatial analysis including buffer, intersection, and searching within a shapefile's attribute table (finding count of a specific attribute occurrence and a single row can have multiple attributes in the same column). This code snippet is very slow. I need to traverse around 2000 records and search for around 179 attributes (count) against 10 different radii.
It takes around 2.45 minutes to traverse against 1 point for 10 radii and search for the count of only 10 attribute occurrences. Is there any way to speed up this process?
I am attaching the code below.
# importing libraries
from openpyxl import load_workbook
from shapely.geometry import Point
import pandas as pd
import geopandas as gpd
#import matplotlib.pyplot as plt
import csv

#loading input excel
book = load_workbook(r'File path of input file.xlsx')
sheet = book.active
bookOut = r'outputfile.csv'
#outputSheet = bookOut.active
#searching for a specific category and returning its count
def bSearch(lst, cat):
     c = 0
     for index, row in lst[0:len(lst)].iterrows():
        if cat in row['Category_t']:
            c = c+1
     return (c)

# creates buffer
def createbuf(p,r):
    bufp = p['geometry'].buffer(distance = r)
    buf = gpd.GeoDataFrame(geometry = bufp, crs = ucs.crs)
    return (buf)

# calculates intersection
def calint(b,u):
    areaoi = gpd.overlay(u, b, how= "intersection")
    return (areaoi)

# shapefiles Input
ucs = gpd.read_file(r'Boundary.shp')
ucs = ucs.to_crs(epsg=32643)

#Business List
business = gpd.read_file(r'List of Categories.shp')
business = business.to_crs(epsg=32643)
businessArray = [
'Cat1','Cat2','Cat3','Cat4','Cat5','Cat6','Cat7','Cat8','Cat9','Cat10'
]
radiusVals=[0.5, 1, 1.5, 2, 2.5, 3.000, 4.000, 5.000, 7.000, 10.000]
# lat long
latlng = []
i = 0
otp = []
for row in range(143, sheet.max_row+1):
    detStore = []
    ind = 0
    for column in "AB":
        cell_name = "{}{}".format(column, row)
        latlng.append(sheet[cell_name].value)
        detStore.append(sheet[cell_name].value)
    for column in "CD":
        cell_name = "{}{}".format(column, row)
        detStore.append(sheet[cell_name].value)
    # PCS SRID
    p1 = Point((latlng[i],latlng[i+1]))
    df = pd.DataFrame({'a':[latlng[i+1],latlng[i]]})
    po = gpd.GeoDataFrame(geometry = [p1], crs = ucs.crs)
    po['geometry'] = po['geometry'].to_crs(epsg = 32643)
    print(i,  "--", latlng[i+1], latlng[i])
    i = i+2
    # calling functions
    for r in range (0, len(radiusVals)):
        outputStore = []
        for d in range (0, len(detStore)):
            outputStore.append(detStore[d])
        outputStore.append(radiusVals[r])
        bufo = createbuf(po,(radiusVals[r]/111))
        aoi = calint(bufo ,ucs)
        busList = calint(bufo, business)
        for b in range (0, len(businessArray)):
            buSearch = bSearch (busList,businessArray[b])
            outputStore.append(buSearch)
        otp.append(outputStore)
print (otp)
fields = ['Location X', 'Location Y', 'Type', 'Name', 'Radius', 'Cat1','Cat2','Cat3','Cat4','Cat5','Cat6','Cat7','Cat8','Cat9','Cat10']
with open(bookOut, 'w',  newline='') as f:
    write = csv.writer (f)
    write.writerow(fields)
    write.writerows(otp)


texts: Problem opening GRASS GUI: libgeos-3.3.6.so: cannot open shared object file:

texts:
I have been using GRASS to process some Landsat imagery but very recently, and as far as I am aware - after having done nothing different or changing the settings, I am trying to use GRASS and on start-up I receive the message:

GRASS 6.4.2 (newLocation):~ > 3D view mode: libgeos-3.3.6.so: cannot open shared  object file: No such file or directory
3D view mode: libgeos-3.3.6.so: cannot open shared object file: No such file or directory
3D view mode: libgeos-3.3.6.so: cannot open shared object file: No such file or directory
g.list: error while loading shared libraries: libgeos-3.3.6.so: cannot open shared object file: No such file or directory
g.list: error while loading shared libraries: libgeos-3.3.6.so: cannot open shared object file: No such file or directory
GRASS module 'g.region' not found. Unable to start map display window.

Consequently the GRASS GUI does not open once I have chosen my location (I have tested a range of locations that previously opened fine).  I am using NX client to run GRASS.

Would anyone be able to help with this error message?


texts:
I have created a python function within field calculator to calculate the time difference between a timestamp (yyyy-MM-dd HH:mm:ss) in a row and the row after (the interval field in seconds). I have used a update cursor with field calculator to do this, however I am receiving null values for both the first and last row.

The first is correctly missed due to the function being reliant on having a value set by the previous row.

I am unsure as to why the last row is being skipped?

The first row of the attribute is formatted as:

The last:

I have had this issue on several other update cursor functions including those on simple single row condition statements.

Am I missing a statement to close the update cursor loop?

My code is:

import arcpy, time, datetime
from time import strftime
from datetime import timedelta, datetime
from arcpy import da

def FindTime(table,date,interval):

  firstRow = True

  with arcpy.da.UpdateCursor(table, ["interval", "date"]) as cursor:

    for row in cursor:
        gap2 = row[1]

        if firstRow == True:
            gap1 = gap2
            firstRow = False
            continue
        timedelta = gap2 - gap1
        row[0] = timedelta.days * 24 * 3600 + timedelta.seconds
        gap1 = gap2
        cursor.updateRow(row)

*Updated from below comment's solution:

import arcpy
def FindTime(fc, datefield, daydiff_field):

    all_dates = [i[0] for i in     arcpy.da.SearchCursor(fc,datefield)]

    diff = [(d1-d0).seconds for d0,d1 in zip(all_dates, all_dates[1:])]
    givediff = iter(diff)

    with arcpy.da.UpdateCursor(fc, daydiff_field) as cursor:
        next(cursor) #uncomment this line if you want to calculate row2 diff as row2-row1, =first row no diff
        for row in cursor:
            try:
                row[0] = next(givediff) #Fetch diffs until list is empty...
                cursor.updateRow(row)
            except StopIteration: #...then break the cursor
                break


texts:
I have this dictionary:

dict1 = ({'ROADS' : ['STREETNAME', 'LEFTFROM', 'LEFTTO', 'RIGHTFROM', 'RIGHTTO'],
            'ITA' :  ['I.ID_BOUND', 'I.BOUND_IDQ']})

I want for every file(ROADS,ITA) to concatenate its columns to one column.

Example:

For the first file concatenate: 'STREETNAME', 'LEFTFROM', 'LEFTTO', 'RIGHTFROM', 'RIGHTTO' to one column in a loop so it can happen for all the other files in the dictionary with their related columns.

I have done this in geopandas but I am having trouble with doing this to Arcpy.

In geopandas:

for k, v in dict2.items():
    df = pd.DataFrame.from_records(data=arcpy.da.SearchCursor(k, v))
    df.columns = v                 # as v are the columns in the dictionary
    df['new'] = df.astype(str).apply(' '.join, axis=1) #new field with the concat columns

I would like this in ArcPy only.

I tried:

In ArcPy:

for k,v in dict1.items():
    with arcpy.da.UpdateCursor(k, v) as cur:
        for row in cur:
            row[2] = row[0] + row[1] #<-should work for all number of columns a file has
            ucur.updateRow(row)

To see how the result should look like:

Example:
if the file has 3 columns a,b,c it will be

from this:

a  b c
1 js 1
1 hd 1
1 jd 1

to this:

a  b c  new
1 js 1  1js1
1 hd 1  1hd1
1 jd 1  1jd1

Update
From BERA's answer:

for k,v in dict2.items():
    arcpy.AddField_management(k,'concat','LONG')
    field_to_update = ['concat']
    with arcpy.da.UpdateCursor(k,v+field_to_update) as cursor:
        for row in cursor:
            row[-1] = ''.join([str(i) for i in row[:-1] if i is not None])

It makes all the rows NULL. Probably a problem with the concatenation of numbers in the columns?
Also is 'LONG' as type correct for the new field?


texts: Is the accuracy or PDOP included in smartphone or GPS camera?

texts:
I would like to properly geolocate picture Exifs on a smartphone or GPS camera and would like to record the accuracy of the GPS position as well.

I have downloaded a few photos from camera review sites and with this short survey I could only find one camera listing the number of satellites.

Any ideas?

Here are the geolocation attributes from my low-end smartphone:

GPSVersionID    2.2.0.0
GPSLatitudeRef  South
GPSLongitudeRef East
GPSAltitudeRef  Above Sea Level
GPSTimeStamp    12:39:37
GPSProcessingMethod GPS
GPSDateStamp    2019:02:01

GPSAltitude 180.4 m Above Sea Level
GPSDateTime 2019:02:01 12:39:37Z
GPSLatitude 32 deg 52' 56.39" S
GPSLongitude    18 deg 38' 32.81" E
GPSPosition (-32.882330556,18.642725000)
            32 deg 52' 56.39" S, 18 deg 38' 33.81" E


texts: Why does my source and target repeats itself when I use pgr_createtopology()?

texts:
I created my path in QGIS v2.18.18, making the lines and a combination of cleaning tool and explode lines tool to get my desired output.
I tried using ROAD GRAPH plugin in order to check whether my path works when getting the shortest path, and it does.
So I used GDAL to put the shp file to my pgadmin, and then using pgr_createTopology() doesnt give me the correct output.
The following is an image of my table.

and here is the the command I ran for it.

SELECT pgr_createTopology('edges', 0.001, 'the_geom', 'gid');


texts: How to create a WCF SOAP web service to find Distance and Direction between two addresses using ArcGIS Online services

texts:
I need to create a WCF SOAP service to find out Distance and Direction between two given addresses. I have to use ArcGIS online service for the same.

Can you please help me out by providing any sample.


texts: Nothing but NAN values using raster interpolation in QGIS

texts:
I am trying to perform a raster interpolation of a contour shapefile with the intent of eventually creating a DEM. No matter what I try, I get NAN values, and nothing is working to fix it. I can't even begin to understand why this is happening, I would love to provide more detail if I had even the slightest idea what the problem was... I've tried changing the cell size, trying points, structure lines and break lines as the type of file, I've tried doing it as a TIN and IDW, changing the coordinate system to UTM (as suggested in another post)... Every possible solution that I can think of.

I've also tried going the grid interpolation route instead of the typical raster interpolation route, no dice. I get a different error with that route though...

I'm clutching at straws, but my guess is that something is wrong with my shapefile, though I can't figure it out. Any suggestions would be appropriated, though I am pretty incompetent with Python if that is the solution...


texts: How do I batch edit field values in my attribute table?

texts: I'm using QGIS.  I have an attribute table that has some very long attribute names because of a join I did, and it seems to be causing errors when I try and select values using an expression and export them as a new shapefile (see the screenshot).  What I would like to do is remove the "CensusACSBlockGroup_" from each column title, but there are a lot of fields. Is there a way to edit the columns and remove this from all of them at the same time?


texts: What do I need to do to add a new QGIS project to qgis-web-client?

texts:
for weeks I'm trying to get QGIS Server with qgis-web-client to work but I'm totally stuck. System is Debian Squeeze. In principle, I get QGIS Server and the qgis-web-client up and running: it's possible to have a look at the example projects "Hello World" and "Natural Earth"; they're loaded in qgis-web-client just fine. The "only" problem I can't get rid of is the fact, that I don't know what to do to have a look at my own projects.

One thing to mention: Actually, we're running our services with NGINX. So I set it up as a proxy for Apache, which ist listening on port 8008.

I've been looking on the internet for a step by step procedure which explains whicht settings I have to set up that my own projects are also working in the web client. As you can see on this website the client is able to render the map and also an extra layer I took from http://www.diva-gis.org .

But saving another project in QGIS (named indonesia_1.qgs) including exactly the same layer and trying to show it in qgis-web-client leads to an empty map.

Any help is really appreciated. I've read the README files and tutorials in the Wiki but I cannot solve this issue alone.

Thanks a lot for your help. I guess I'm just doing a tiny mistake here -- somewhere.

Niklaas


texts:
I have set of polygons whom I would like to check their compactness.
I'm working in Jupyter Notebook and my geometry is GeoPandas.
I would like to try Polsby-Popper test and Schwartzberg test (and other test that are here:  https://fisherzachary.github.io/public/r-output.html ) but I can't find any simple way to do it.
I don't find way to create the necessary circles, e.g "circle whose area is equal to the area of the polygon " or "circle whose circumference is equal to the perimeter of the polygon ".
I have found script to do find smallest enclosing circle but it seems like it works with points and I have polygon in GeoPandas  (https://www.nayuki.io/res/smallest-enclosing-circle/smallestenclosingcircle.py).
If anyone knows any library/package/ any idea of how I can create circles from given perimeter or area. It's important to say that I have around 70k polygons to check so I also have the memory issue.
My end goal is to get ideas how can I calculate this in Python in Jupyter Notebook.
Edit: I have found this library but struggling with install it/use it
(https://jblindsay.github.io/wbt_book/available_tools/gis_analysis_patch_shape_tools.html)
(https://github.com/jblindsay/whitebox-tools/blob/master/src/tools/gis_analysis/related_circumscribing_circle.rs)


texts: Colour Scale Bar for One-Band Rasters Error

texts:
I have a vector layer of a map with a number of regions. I have coloured those regions using a gradient blue colour. I would like to create a colour scale bar for this map. I downloaded the Colour Scale Bar for One-Band Raster plugin from this website: http://www.bc-consult.com/free/bccscbar.html.

I created a .qml file by going into Styles and clicking on Save Styles - QGIS Layer Style File, but when I try to open this file in the plugin, it gives me an error saying "Error Wrong Input File (must be a colour table .txt or a custom QML file!).

Maybe the error is to do with the fact that I made the .qml file using a vector map and not a raster map. Could someone please help me in this. Is there a way to make the colour table in .txt from Styles (in a vector map) and then opening that file in this plugin to make a Colour Scale bar.


texts: Free printable background map: land use, rivers, etc

texts:
I am looking for a free printable background map of europe which should look more or less like the OCM landscape:
land use, major rivers, coast lines
Is there any WMS or raster set whiche enables me to make quick distributions maps of species in Qgis and print it.

Borders would be useful but not essential, labels of cities, major roads etc. are not need nor wished.


texts: Updating geometries using another table whose records need to merge in SQLDeveloper

texts:
I'm looking for some help with the following.
I have 3 tables in my Oracle Database.
Table A:
IDnumber  Geometry
001       Geom1
002       Geom2
003       Geom3

Table B:
IDnumber  Fieldnumber (neither of these are unique fields)
001       Field1
001       Field2
002       Field3
003       Field1
003       Field4

Table C:
Fieldnumber  Geometry
Field1       Geometry1
Field2       Geometry2
Field3       Geometry3
Field4       Geometry4

I want to update the geometries from Table A with the geometries from Table C using Table B as the connector. If there is more than 1 Fieldnumber for each IDnumber in Table B, I want the geometries from Table C to merge.
My update SQL so far:
Update TableA a
Set
a.Geometry = (select c.Geometry from TableC c where c.Fieldnumber = b.Fieldnumber)
Where
a.IDnumber in (select b.IDnumber from TableB b)

This only works when the TableB IDnumbers are unique. How can I merge them? With StUnion? If so how does that work? I don't understand the code.


texts: Windows software to convert GeoJSON (> 6Mb) to Excel/CSV?

texts:
I am trying to convert OpenStreetmap data GeoJSON to Excel (or .csv files). Few online converters are there.

However, when the file is big (example: 6Mb Geojson), the browser is becoming unresponsive. So I am looking for converters that will handle big files. My current software, Opal converter, is unresponsive when converting big files.

Is there any windows desktop software available to convert GeoJSON greater than 6Mb to Excel/CSV?


texts: Trying to find percentage of overlap using turf intersect returns 'found non-noded intersection between LINESTRING'

texts:
This worked just the other day and now it's throwing an error for every call.

The original code (before adding the try/catch) is this:

var cbgShape = turf.polygon(shapeData.geometry.coordinates.geometry.coordinates);
var intersection = turf.intersect(isochronePolygon, cbgShape);

Then, it started throwing the 'non-noded intersection' error from turf.intersect, so I wrapped it in the try/catch.

//determine intersection between the current polygon and the isochrone polygon
var cbgShape = turf.polygon(shapeData.geometry.coordinates.geometry.coordinates);
var intersection;
try {
    intersection = turf.intersect(isochronePolygon, cbgShape);
} catch (e) {
    console.log(e);
    intersection = null;
}

var cbgArea = turf.area(cbgShape);
var overlapArea = 0;
var overlapPercentage = 0.0;
if (intersection) {
    overlapArea = turf.area(intersection);
    overlapPercentage = (overlapArea / cbgArea) * 100;
}

And now, there are no intersections at all, which is of course, incorrect.  Many of the cbgShape polygons are completely within the isochronePolygon.

This all worked just fine a few days ago and the code that generates the isochronePolygon (which is from Mapbox) hasn't changed, and the code that generates each of the cbgShape polygons hasn't changed.  The data that is used for the cbgShape polygons hasn't changed (it comes from static data stored in the database).

So, where do I begin look to resolve this?  I have already tried reverting to older code, but it's still happening.  Help!


texts: Which algorithm (Vincenty or Haversine) is used to calculate distance in GeoTools?

texts:
I want to use GeoTools library to calculate geographic distance by using GeodeticCalculator (as mentioned in here http://docs.geotools.org/stable/userguide/library/referencing/calculator.html ).

Does this library use Vincenty formula, Haversine, or something else?


texts: Using ArcGIS REST Service to query particular features inside BBOX?

texts:
I really have 0,0000... knowledge about OGC web services (WMS / WFS ; only limited basics). Therefore maybe this 'stupid' question.

My question might be misplaced / or have the wrong tags... But it's hard to determine what is what without knowledge of it all...

I have a WMS - I use it for retrieving a map image according a specific BBOX.

  https://geoservices.informatievlaanderen.be/raadpleegdiensten/ANB/wms?REQUEST=GetCapabilities&service=WMS

Since they don't have a matching WFS I'm not able to get the features inside the BBOX. Since the WMS only deals with a X/Y point for getting a feature on that particular point. I have tried with the BUFFER parameter but does not seem to work.

I also tried rendering the whole BBOX to a 1/1 (width/height) image and getting the features for the single X/Y point...  But those results don't always seem correct...

So I gave up on that one, until I found following ArcGIS REST Service...
Is there a possibility to obtain the features inside a specific polygon or bbox using the link below? Opening the link gives me the impression there is a query possibility...

  https://www.geoloket.be/gwserver/rest/services/Milieu/AGS_historische_bossen/MapServer/2/query


texts:
I'm using the Field Calculator .replace function to replace WEST, EAST, NORTH, SOUTH to W, E, N and S and I'm using 4 scripts !PREDIR!.replace("WEST*","W") and so forth.

I would like to use one string to do this, how can I do this?


texts: How can i an add image overlay depending on the options in the polygon properties in a leaflet map?

texts:
I have a multi-polygon geojson file ( mapInfo ) that has one of the properties named "field_hazards". The hazards include "gale/storm, heavy rain, thunderstorm and Freezing rain "How can i an add image overlay depending on the field_hazard options?
Note the geojson file is got dynamically as a drupal feed.
I have made the following code but it returns an empty map, the polygons are not rendered.

function getImage(d) {
  return d === 'Freezing Rain' ? "http://mymap:8082/images/weather-images/43n.png" :
    d === 'Thunderstorm' ? "http://mymap:8082/images/weather-images/11.png" :
    d === 'Heavy Rain' ? "http://my:8082/weather-images/02n.png" :
    d === 'Gale\/Storm' ? "http://mymap:8082/images/weather-images/15.png" :
    "http://my:8082/images/weather-images/09.png";
}
var imageUrl = getImage(feature.properties.field_hazards);
var imageLayer = L.imageOverlay(imageUrl, imageBounds).addTo(map).bringToBack();
var boxOptions = {
  fillOpacity: 0,
  opacity: 0,
  onEachFeature: onEachBox
};
//create the image interaction box
var imageBox = L.geoJson(mapInfo, boxOptions).addTo(map);
//zoom in to fit GeoJSON layer
map.fitBounds(imageBox.getBounds());


texts:
I am carrying out a piece of work in ArcGIS 10.4.1 that involves realigning polygons to an updated vector map. The polygons in question are planning permission boundaries that are created showing the outline of the property the planning permission request refers to. Over a number of years a single property can have had multiple polygons drawn to reflect the different planning applications made against it. These polygons have been drawn by various different people over the years and depending on how accurate they were whilst creating them, they may have missing or extra vertices or no longer match the updated vector map due to changes in the the most recent mapping updates.

I am looking for a quicker way to achieve this than editing the individual polygons. Ideally what I would like to do is define a polygon shape created from my new base map layer and replace the geometry of all the planning permissions in one go. I do not know if this is possible but I have approx 60,000 polygons that potentially need replacing / realigning to new geometry and I cannot see an end in sight doing these individually!

Edited to include screenshot, In this screenshot, the red polygons is the planning layer, I need to align these polygons with the geometry of the green vector map, in this screenshot I have four misdrawn polygons that need to be redrawn to match the green vector map....


texts: Getting ERROR 010161 in geopackage using ArcMap

texts:
I am using ArcMap 10.7.1.
I'm trying to share a model via geopackage. The model works fine while open in my .mxd file and shares without any errors but when I open the newly saved geopackage I receive error 010161 once it gets to the reclassify tool. I can't figure out why it only fails here. Seems to be an issue with the pathway but I don't understand why the previous tools work if that is the case.
Also, the geopackage should have bundled all the data so file navigation should only come from the scratch workspace?
Also the same happens for each of the reclassify steps.


texts: Choosing the right projection using PyProj

texts:
I'm starting with GIS work trying to analyze movement of cars on a track.
I've understood that in order to use a library like shapely I need to transform from GPS spherical coordinates to Cartesian and for that I need a projection.
I'm confused with how to choose the right projection. The area of interest this way or another is not big but I couldn't wrap my hand around which projection to choose from the "catalog".


texts: How do I make a ratio image on Erdas Imagine 2014?

texts: I have a stacked image but want to divide bands 4/5/. I have built a model on the model maker but the resulting image is completely black after running the model. When I open the attributes table the two categories are both "1" and are black


texts: Is there any way to reset fid value of data table? I'd to remove all data from table and I'd like to got '1' in fid column in first row for new entry.


texts: Calculating height in raster layer pixel wise using QGIS C++ API

texts:
I am making a custom QT application using QGIS 3.14 C++ APIs.
I am able to render a raster input file(.Tiff). But sometimes Raster file can be of DTED type(with .tiff extension)(DTED file).
How can I calculate elevation data pixelwise and display it on QgsmapCanvas?
Which API should be called while hovering (or by clicking somewhere on the layer) on QgsmapCanvas.
I need suggestions from the developer's perspective?
I am not using QGIS standalone application.


texts:
How to create a map with polygons from all the different layers?

I'm using ArcGIS Pro and it does not have Overlay analysis tool.


texts: Setting anchor/pivot point to layer features (point and text) by which symbol or text will rotate in ArcMap?

texts:
I know I can change the rotation and pivot point when editing, but if I have point layer that contains for example an arrow symbol (EMF picture) and I need the arrow to rotate around the bottom part of the arrow, can I somehow set it?

In QGIS for a text layer, there is an option Quadrant (Properties - Labels - Placement) where you can choose values from 0 - 9. Each number represents a position (pivot point). And for points there is an option anchor point (Properties - Style) and there are 5 options for the anchor (pivot) point.

Is there anything like that in ArcMap 10.2.2 ?

This is how I did it in QGIS. That arrow is the EMF file and I've set the anchor point to center and right so if there is a rotation, the arrow will be rotating around the right part and that's what I need to do in ArcMap. And I also need to do that to many different symbols and texts.


texts: Dissolve polygons - group to others with touching corners

texts: I am working with some data that shows land management.  There are many places where corners "touch".  I would like to dissolve my dataset so that features are grouped (dissolved) with those that they share a touching corner.  Each section (square) is a unique feature (row).  In the image below I've circled the features that I would like to group in this dissolve.  I am using ArcGIS and have an editor and SA license.


texts: I've made a regular grid and regular points for a grid, I want to set this inside of a boundary so I have deleted the points outside the boundary. However, when I do this the ID still shows up counting the deleted point. I want the visible points to start with and ID of 1 or 0.


texts: Append new data to existing feature class using Load Data in ArcGIS Pro

texts:
I am trying to append new rows of data for a feature class that shares the exact same fields (column headers and types) but wasn't able to using Load Data in ArcGIS Pro. In ArcMap, you could easily hit reload table cache and the table will update by calling the CSV that it uses to plot the feature class. But it seems very difficult to do it in ArcGIS Pro.
When using "Load Data" function, I tried 2 options:

Input schema must match target schema - the input table has the exact same header with new data to be appended. It kicked Error 000466.
Use the Field Map to reconcile schema differences - Same input table, but I also double checked in the Field Map section to make sure the properties of each output field is the same as what is in the feature class now. This step gave the Error 001156.

I might not be exactly sure how to use this tool.


texts:
I am a little stumped with determining the area of polygons in the python console. I want the output to be in meters^2, so using:

qgis.utils.iface.mapCanvas().setMapUnits(0)

qgis.utils.iface.mapCanvas().refresh()

cLayer = qgis.utils.iface.mapCanvas().currentLayer()

for f in cLayer.getFeatures(): print f["ID"], f.geometry().area()

I get the the output of:

1 0.000114895410903
2 0.00014639997019...  - which suggests to me that the map units are still in in decimal degrees and not meters squared like I expected.

with some background reading, I then specified a localised CRS:

wkt = 'PROJCS["SACRS_NO_27",GEOGCS["WGS 84",DATUM["unknown",SPHEROID["WGS84",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",27],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
 and applied it with:

coo = QgsCoordinateReferenceSystem(wkt) and
cLayer.setCrs(coo)

I can check the mapUnits of cLayer with
print cLayer.crs().mapUnits() and this gives me '0' for metres.

Yet for f in cLayer.getFeatures(): print f["ID"], f.geometry().area() still gives me the same output as before, as if the units are still in decimal degrees.

What am I missing here?


texts: How install postgis_raster extention on PostgreSQL 11.5?

texts:
I want install postgis_raster extension on PostreSQL 11.5 + PostGIS 2.5.3 + MacOS 10.14.6

There is a install script on https://postgis.net/install/.

All installations of extensions work fine, except postgis_raster:

CREATE EXTENSION postgis_raster;
2019-11-12 17:39:19.416 CET [88243] ERROR:  could not open extension control file "/usr/local/share/postgresql/extension/postgis_raster.control": No such file or directory

What can I do?


texts: How to store Google Maps DrawingManager Rectangle overlay in PostgreSQL Postgis GEOMETRY(GEOMETRY, 4326)

texts: What is the correct way to store a Google Maps DrawingManager rectangle overlay? The API basically provides me with NE & SW lat/long coordinates. Should I use BOX()? Is it possible to create a BOX() from a GeoJSON object?


texts: Add "Percent Complete" chart based on two feature classes to Operations Dashboard

texts:
I am trying to create an Operations Dashboard in ArcGIS Online and want to track progress towards a goal.  The goal is a sum of features (count of potential_customers, subtyped as residential).  Progress is a sum of counts from a separate feature class of building_type where subtypes "single family" count once, but "duplex" count twice.

I would like to display a number representing percent complete, or maybe a graph depicting progress next to goal. How do I accomplish this?


texts: I have a point dataset which I'd like to Krige, ideally using an open-source software package. If possible, I'd also like to choose the semi-variogram model during the process to improve the estimation.


texts:
I've been trying to convert all my ArcGis projects into QGIS 3.4.10.
I have a problem with my raster data sets. They're all .tiff data and have the right projection but they still won't fit to my shapefiles and wms servers.

Do I need the .tfw data or something else?


texts:
I am using QGIS 2.8.1-Wein 32bit version.  When running QGIS, it does not open. After displaying like below the picture, nothing happens. If you can solve this problem, please give me a solution.

I also tried re-installation.


texts: Create ESRI shapefile with attributes from mySQL

texts:
I have some geometry which can be perceived as regions within a country and wish to create a shapefile which contains all the regions along with some attributes (e.g. name, population, area code, etc.)

I have the geometry from the regions in x,y WGS84 lat lon (not polygons so cannot be sure they are closed loops which is priority) currently stored in mySQL.

I wish to be able to create a single shp (with associated dbf and shx files) to contain this data.

Any pointers on command line tools and methods to do this would be gratefully received.


texts:
I am trying to create an animation using Time Manager for multiple vector layers.

Each layer is a graduated map of rent prices for a quarter of a year.

I have 8 quarters, each in a separate layer. 4 for 2016, 4 for 2017. I would like to run an animation that shows rent prices in Q1 2016 and runs until Q4 2017.

My layers have no date column. So I imported a csv with dates for each quarter, e.g. for Q1 2017 there is a column Q1_2017 with the value 2017-01-01 and I computed a join to this csv for each of the 8 vector layers to give them a unique date column.

When I run time manager I am adding a layer with the following settings:

When I add all my layers in a similar way and try to run the animation I get an error like the below:

Is it my settings or is it the join?

Does the date field need to be in the vector layer itself?


texts: Is it possible to change automatically the font size of the labels which are not shown in the current map canvas on print composer?

texts:
I would like to make a rule-base labelling depend on: is a label shown in the current view or not. ("Show all labels for this layer (including colliding labels)" is unchecked.)

I would like to use it on a polyline layer, with curved placement.

For example:
The default font is 8p. If a label is not shown on the current view (because of colliding labels, I think), then the font size of that label changes to 7p, so it may come to be shown. If it does't, the font size of it changes to 6p, etc.


texts: What's up with the GeoJSON spec and CRS as a IRM?

texts:
I'm writing a web API to our GIS data, currently revamping our geocoding service. I'm creating the option for folks to specify the output type as GeoJSON so I'm trying to follow the spec. Our data is stored in nad83 datum UTM zone 12 north projection. The GeoJSON spec says to add a CRS. The wkid is 26912 but what is the OGC CRS URN?

I can understand that the OGC CRS URN is broken down into a few parts.

urn:ogc:def:crs:OGC:1.3:CRS84

urn is the identifier, ogc is the organization, def is another static deal, crs is the type (coordinate reference system), OGC is the authority, 1.3 is the version and CRS84 is the projection.

Why is the authority of utm 26912 the European Petroleum Survey Group?

Do any mapping frameworks even use this CRS name?

I want to follow the spec but it seems that more mapping frameworks are using the wkid.

http://spatialreference.org/ref/epsg/26912/

seems like the best place to get the info but they don't even list the urn. Is the GeoJSON spec just silly or what?


texts: How to filter two map layers with one button on CartoDB?

texts:
I am working on a mapping project in which I want the end user to be able to filter the data. This tutorial (http://docs.cartodb.com/tutorials/toggle_map_view.html) gave me a way of doing this, such that the end user can effectively filter by clicking on buttons. Where I am now stuck is that my map has two layers (one for points, one for lines), and I want the button to apply to both layers rather than only one. Which parts of the code do I need to change in order to make this happen, and what changes do I need to make?

The part that I am looking at in particular from the tutorial is:

var sublayers = [];

  cartodb.createLayer(map, layerUrl)
  .addTo(map)
  .on('done', function(layer) {
    // change the query for the first layer
    var subLayerOptions = {
      sql: "SELECT * FROM ne_10m_populated_places_simple",
      cartocss: "#ne_10m_populated_places_simple{marker-fill: #F84F40; marker-width: 8; marker-line-color: white; marker-line-width: 2; marker-clip: false; marker-allow-overlap: true;}"
    }

    var sublayer = layer.getSubLayer(0);

    sublayer.set(subLayerOptions);

    sublayers.push(sublayer);
  }).on('error', function() {
    //log the error
  });

  //we define the queries that will be performed when we click on the buttons, by modifying the SQL of our layer
  var LayerActions = {
    all: function(){
      sublayers[0].setSQL("SELECT * FROM ne_10m_populated_places_simple");
      return true;
    },
    capitals: function(){
      sublayers[0].setSQL("SELECT * FROM ne_10m_populated_places_simple WHERE featurecla = 'Admin-0 capital'");
      return true;
    },
    megacities: function() {
      sublayers[0].set({
        sql: "SELECT * FROM ne_10m_populated_places_simple WHERE megacity = 1",
        //as it is said, you can also add some CartoCSS code to make your points look like you want for the different queries
       // cartocss: "#ne_10m_populated_places_simple{ marker-fill: black; }"
      });
      return true;
    }
  }

Edit: Right now it seems that I could accomplish this by just making the button functions do two actions (i.e. select x from one table, select x from the other table). Is this possible? If so, how?


texts:
Based on the previous questions I asked, I got a way to fetch the clicked point information, but I want the exact value of gray_index, but as you can see the answer is text or HTML. In the image below, you can see the answer in the console as below.

Gere is my layer code :
const UnTiled = new ol.layer.Image({
        opacity:0.9,
        title:'UnTiled',
        source: new ol.source.ImageWMS({
            ratio: 1,
            url: 'http://localhost:8080/geoserver/world/wms',
            params: {
                     'VERSION': '1.1.1',
                  "STYLES": '',
                  "LAYERS": '       world:world',
                  "exceptions": 'application/vnd.ogc.se_inimage',
            }

          })

        });

and here is my function for fetching data :
map.on('singleclick', function(evt) {
    document.getElementById('info').innerHTML = '';
    var view = map.getView();
    var viewResolution = view.getResolution();
    var url = UnTiled.getSource().getFeatureInfoUrl(
        evt['coordinate'],
        viewResolution,
        'EPSG:3857', 'EPSG:4326',
        {'INFO_FORMAT': 'application/json'}
      );
      console.log(url);
      if (url) {
        fetch(url)
          .then(function (response) { return response.text(); })
          .then(function (html) {
            html;
            console.log(html)
          });
      }

     });

The URLs are like:
http://localhost:8080/geoserver/world/wms?0=E&1=P&2=S&3=G&4=%3A&5=4&6=3&7=2&8=6&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=%09%09world%3Aworld&STYLES=&LAYERS=%09%09world%3Aworld&exceptions=application%2Fvnd.ogc.se_inimage&X=50&Y=50&SRS=EPSG%3A3857&WIDTH=101&HEIGHT=101&BBOX=5203235.632400199%2C3561089.2605374325%2C5697324.583235578%2C4055178.2113728123


texts: QGIS Display features by attribute value

texts:
I have a layer called "streets".  It has an attribute called "traveled".  The values are either "yes" or "no".  If the value is "no", I would like to display those streets in a different color (than the ones with "yes" as the value).

Also this "streets" layer has an attribute called "miles".  If the "traveled" attribute is "yes", I would like to retrieve the total "miles" from the "traveled" streets only.


texts: Show popup on marker click on static image in Openlayers

texts:
Im complete newbie in scripting and I have problem which I can't solve (and I searched whole web for days for coming here for answer). Im trying to make very simple completely static image map. I managed to put markers on it and the problem I encountered is that, even I used lot of examples from the web, I can't make popup to open on click on marker which will then contain some description for that specific marker.

For the simplicity, Ill put the usual 'Online communities' static image map with one marker in the code.

index.html

<!DOCTYPE html>
<html>
<head>
    <title>Example</title>
    <script src="https://code.jquery.com/jquery-1.11.2.min.js"></script>
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com   /bootstrap/3.3.4/css/bootstrap.min.css">
    <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/ol3/3.6.0/ol.css" type="text/css">
    <script src="https://cdnjs.cloudflare.com/ajax/libs/ol3/3.6.0/ol.js">    </script>

</head>
<body>
    <div id="map" class="map">
    <script src="index.js"></script>

</body>
</html>

index.js

import 'ol/ol.css';
import Map from 'ol/Map';
import View from 'ol/View';
import {getCenter} from 'ol/extent';
import ImageLayer from 'ol/layer/Image';
import Projection from 'ol/proj/Projection';
import Static from 'ol/source/ImageStatic';

  // Map views always need a projection.  Here we just want to map     image
    // coordinates directly to map coordinates, so we create a projection that uses
    // the image extent in pixels.
    var extent = [0, 0, 1887, 949];
    var projection = new ol.proj.Projection({
        code: 'xkcd-image',
        units: 'pixels',
        extent: extent
    });

    var iconFeature = new ol.Feature({
    geometry: new ol.geom.Point([586, 256]),
    name: 'Null Island',
     population: 4000,
    rainfall: 500

});

var iconStyle = new ol.style.Style({
    image: new ol.style.Icon({
        src: 'https://i.ibb.co/z6HKFRB/iconfinder-map-marker-299087-1.png'
    })
});

iconFeature.setStyle(iconStyle);

var vectorSource = new ol.source.Vector({
    features: [iconFeature]
});

var vectorLayer = new ol.layer.Vector({
    source: vectorSource
});

var map = new ol.Map({
layers: [
    new ol.layer.Image({
        source: new ol.source.ImageStatic({
            attributions: [
                new ol.Attribution({
                    html: 'Online communities - example - xkcd'
                })
            ],
            url: 'https://imgs.xkcd.com/comics/online_communities.png',
            projection: projection,
            imageExtent: extent
        })
    }),
    vectorLayer,
],
target: document.getElementById('map'),
view: new ol.View({
    projection: projection,
    center: [500, 500, 500, 500],
    extent: [500, 500, 500, 500],
    zoom: 2.7,
    minZoom: 2.7,
    maxZoom: 2.7
})
});


texts: How to delete fields from a Geopackage attribute table?

texts:
Using QGIS 2.14.0 on a Microsoft Surface Book with Windows 10, I created a geopackage from a layer using Save As...

I can successfully add a field to the geopackage's attribute table, calculate changes to its values, and then save the changes.

However, when I try to delete that field, QGIS responds with "Provider does not support deleting fields".  Sure enough, the field cannot be deleted.  Please see attached image.

If I repeat the process, instead using either a shapefile or a spatialite layer, I can delete the attribute field.  So the problem seems to be with the geopackage.

In summary, how do I delete a geopackage attribute table field?


texts: In_memory workspace in geoprocessing services

texts:
I have a arcpy python script that runs a process that I publish as a geoprocessing service on ArcGIS Server 10.3.1.  The python script uses various input feature classes. These feature classes are within geodatabases that are registered on the ArcGIS Server.

My python script also uses the in_memory workspace for storage of some temporary data.  The in_memory feature classes are not the input or final output, there are only intermediate  data.

When trying to publish the geoprocessing service, the publishing tools want to copy the in_memory feature classes to the server, since the in_memory workspace is not registered with ArcGIS Server. The specific warning I get is Data source used by Script MyToolThatDoesStuff  is not registered with the server and will be copied to the server: in_memory\myTempData

I let the publishing tools copy the in_memory data to the server and the my geoprocessing service works as expected.  However, I suspect that every time the tool is run, the in_memory data is copied to the arcgisserver\directories\arcgissystem\arcgisinput\ MyToolThatDoesStuff  .GPServer\extracted\v101\data1.gdb geodatabase and never removed. Over time this gdb bloats, slows down the geoprocessing service and ultimately fills the disk, creating mass problems.

My questions:

Is there a way to prevent the geoprocessing tools from publishing
the in_memory workspace to the server?
Is there a way to register the in_memory workspace with ArcGIS Server?
Is using the in_memory workspace in scripts that are published as geoprocessing services not “best practice”.  If that's the case, how should temporary, intermediate data be processed in arcpy scripts?


texts: How to create graphic margin around labels in GeoServer using SLD?

texts:
I have a little problem using VendorOption - graphic-margin around labels of roads. In some cases they're drawn and some places they're not. You can see it on this picture. In the square its drawn correctly but in the circle the graphic margin is missing.
There is my code under the picture.

Any idea how to solve this?

<Rule>
      <MaxScaleDenominator>10000</MaxScaleDenominator>
      <Name>S1p</Name>
      <Title>silnice 1. tridy paprsek</Title>
      <ogc:Filter>
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>TYPSIL_K</ogc:PropertyName>
          <ogc:Literal>S1p</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
  <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#686868</CssParameter>
          <CssParameter name="stroke-width">4</CssParameter>
          <CssParameter name="stroke-opacity">1</CssParameter>
        </Stroke>
      </LineSymbolizer>
  <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#FF7F7F</CssParameter>
          <CssParameter name="stroke-width">2</CssParameter>
          <CssParameter name="stroke-opacity">1</CssParameter>
        </Stroke>
      </LineSymbolizer>
      <TextSymbolizer>
   <Label>
     <ogc:Literal>I/</ogc:Literal>
     <ogc:PropertyName>SILNICE</ogc:PropertyName>
   </Label>
   <Font>
       <CssParameter name="font-family">http://www.edpp.cz/arial.ttf</CssParameter>
       <CssParameter name="font-size">12</CssParameter>
       <CssParameter name="font-style">normal</CssParameter>
       <CssParameter name="font-weight">normal</CssParameter>
   </Font>
   <Graphic>
      <Mark>
        <WellKnownName>square</WellKnownName>
        <Fill>
      <CssParameter name="fill">#FF7F7F</CssParameter>
        </Fill>
      </Mark>
    </Graphic>
   <LabelPlacement>
     <LinePlacement/>
   </LabelPlacement>
   <Halo>
     <Radius>2</Radius>
     <Fill>
       <CssParameter name="fill">#FF7F7F</CssParameter>
     </Fill>
   </Halo>
   <VendorOption name="followLine">true</VendorOption>
   <VendorOption name="repeat">500</VendorOption>
   <VendorOption name="graphic-resize">stretch</VendorOption>
   <VendorOption name="graphic-margin">2</VendorOption>
 </TextSymbolizer>
</Rule>


texts: Do I need to explicitly release extension licenses?

texts:
I'm writing an application using ArcObjects. I initialize the license with code similar to this:

var extensionCodes = new[] {
    esriLicenseExtensionCode.esriLicenseExtensionCodeNetwork,
    esriLicenseExtensionCode.esriLicenseExtensionCodeSpatialAnalyst
};

lic = new AoInitialize();
lic.Initialize(esriLicenseProductCode.esriLicenseProductCodeAdvanced)
foreach (var ext in extensionCodes) {
    lic.CheckOutExtension(ext);
}

(I have more result checking and error handling in my actual code, as well as a pre-check for availability.)

All the samples and documentation I can find say to use Shutdown to release the licenses:

lic.Shutdown();

But none of the examples I could find used extensions. Is Shutdown enough to release them back to the license server, or do I need to check those back in explicitly like so:

foreach (var ext in extensionCodes) {
    lic.CheckInExtension(ext);
}

I'm using ArcGIS 10.3.1 currently, but I would like to know if previous versions are known to be different.


texts: Unable to install ECW support on lubuntu 14.04

texts:
I try unsuccessfully to install the ECW support (mainly to use ECW file in QGIS 2.2).
We attempt to follow and adapt the howto provided by makina-corpus.
My first supprise was the repository ubuntugis-unstable is unavailable for trusty and the repository http://qgis.org/debian doesn't have libgdal-ecw-src package.
Finally in desperation we add ubuntugis-unstable for the raring distribution in source.list.

And apply the following commands

chmod +x ERDAS-ECW_JPEG_2000_SDK-5.1.1.bin
./ERDAS-ECW_JPEG_2000_SDK-5.1.1.bin # Select option desktop read only
sudo cp -r ERDAS-ECW_JPEG_2000_SDK-5.1.1 /usr/local/
sudo ln -s /usr/local/ERDAS-ECW_JPEG_2000_SDK-5.1.1/Desktop_Read-Only/lib/x86/release/libNCSEcw.so /usr/local/lib/libNCSEcw.so
sudo ldconfig
sudo apt-get install libgdal-ecw-src
sudo gdal-ecw-build /usr/local/ERDAS-ECW_JPEG_2000_SDK-5.1.1/Desktop_Read-Only

The ECW format seems to be installed

gdalinfo --formats | grep -i ecw
  ECW (rw+): ERDAS Compressed Wavelets (SDK 5.1)
  JP2ECW (rw+v): ERDAS JPEG2000 (SDK 5.1)

But .... Gis crash when a ECW raster file is open with the following message :

QGIS died on signal 4Abandon (core dumped)

If anyone have an idea


texts: How can I add r.local.relief to the QGIS Processing Toolbox?

texts:
I'm working with QGIS 2.12 and GRASS 7.0.3RC1. I'm making a model using the model builder and trying to add the GRASS addon "r.local.relief" in my QGIS model, but the addon does not show up in the Processing Toolbox of QGIS. The addon installs sucessfully, and the file for it is shown in the folder with the other GRASS tools as a '.bat' file (the other files are '.exe'). Is there anyway that this tool can be added to the processing toolbox of QGIS 2.12 to be used on the model builder?
I've searched all afternoon but I can't seem to find any pages on the integration of that tool on QGIS, only on the use of it in the GRASS GUI.