texts: Creating a line that goes through two point points using python (used for creating a utility line where there are two pole structures)

texts:
I am starting up an enterprise GIS program at the Utility company I work for.  We are a Generation and Transmission company.  Our data was bad when I got here, it was in NAD27 in a .dwg and was off by at least 10 meters.

After a lot of hard work we got all the poles rectified and imported into a Geodatabase so they match aerial photography so everything lines up in ArcGIS Portal.  Lines had already been created in the same coordinate system, but they too are considerably off.  We can use the point to line geoprocessing tool, but there is a problem.  We have a lot of what are called H-frame structures (they usually are higher voltage transmission lines) and each pole has an entry in our oracle database as they are considered separate structures.  So we will have a pole 0222 and pole 0222A, about 15 feet apart depending on the construction.  If they were a single pole structure I could the point to line tool and everything would be peachy.  But because we have the two pole structures when we run the tool we end up with something that looks like this.

We, of course, don't want this zig zag type of line and want it to run between the two structures so it looks like this

What I'm trying to do then is write a script where it goes through the entire pole dataset (we have 30,000+ poles so drawing them in is out of the question) and selects a pole that is within 50ft (the span between each structure is much longer than that so it will select the other pole in the structure).  I have it set to calculate the UTM coordinates of the structure adds the northings and eastings together and divides by 2 to find the point in the center between them.  I suppose I could also tell it to draw a line in_memory and calculate the midpoint, but my approach seems more simple (maybe it isn't).  Anyway, at this point I'm stuck using a for loop with the search cursor to select the line and just print out the easting and northing of the lines, then I will move on to the math part and the creating lines based on those midpoints.  I had someone here who knows a lot more about python take a look at my code, but he wasn't familiar with the arcpy library so he said he couldn't help much.

This is the point that I'm stuck at.

I'm pretty rusty on my python and am not a programmer by trade.

#This Script will convert a series of poles to lines#
#Power Co-op. Inc. Format"
#Peter#
#12 October 2017"

#Import Modules
import arcpy
from arcpy import env
import os
import math
import numpy

#Set Workspace

arcpy.env.workspace = "in_memory"

#set Parameters#
fc =r"P:\Scripts\pole_2_line\demo_data.gdb\demo_poles_short"
fc_layer = arcpy.MakeFeatureLayer_management (fc, 'demo_poles')
pole = "POLE_NUMBE"
line = "LINE_NUMBE"
coordinates = ['SHAPE@XY']

#select adjacent lines
with arcpy.da.SearchCursor (fc_layer, coordinates) as cursor:
    for row in cursor:
        print row
        near_pole = arcpy.SelectLayerByLocation_management (fc_layer, 'WITHIN_A_DISTANCE',fc_layer,"50 feet", "NEW_SELECTION","NOT_INVERT")
        print near_pole


texts: returning exceptions to user via AGS REST

texts:
I've got a python script that I run as a standalone and as a gp service in AGS (moving to 10.1 now...). For the standalone execution, I'm throwing a custom exception if there is no data found based on the user's input:

    arcpy.MakeTableView_management(lookupTable, luTV, '"COMID"='+gage)
    if arcpy.GetCount_management(luTV)==0:
        raise Exception('UnknownIdentifier')

I'm expecting users of the AGS service to be scripting against the REST API, most likely with python, and I want to be able to pass back the same exception through the AGS web service. I have users that plan to hit my service hundreds of times using some kind of iteration--I need a clear way to let them/their code know that a request had a problem.

I'm having trouble finding information on how to I should modify my script. I can see from this question that there will be some differences if the initial request is executed synchronously (somehow directly passed) or submitted as an asynchronous job (use GPResult), although I didn't understand what these differences mean for either the user's or my code.

I'm not seeing what I think I need here, although I probably (definitely) don't understand it all. Don't know if I have to fiddle with the actual configuration of AGS's REST API--I'm really hoping to avoid needing admin privileges.

Is it really as simple as ditching my raise Exception() syntax and using arcpy.AddMessage() instead?

Thanks for any guidance!

Updates

related posts:

Setting Output Parameter in ArcPy for Server Tool
Using arcpy.SetParameter() with Python Script Tool published as Geoprocessing Service?
Crucial concepts for (ESRI) Geoprocessing Services


texts:
I would like create something like that in QGIS:
http://2.bp.blogspot.com/--EFSaghZ92s/VqQXXIoE4AI/AAAAAAAABKw/KGkoB8zZ6K8/s640/scs2%2Bmap.jpg

I have a huge amount of points and I can import in QGIS (lat, long). I can create also heatmap but the result isn't like image above.
if it is useful to the purpose, I have also the possibility to import a huge amout of line (I can join all my points with index before import in QGIS).

How can I obtain the lines density like the image above?


texts: QGIS 3.14: Remove intersect points on a line

texts:
Is there a way to remove the point that intersect with the red line instead of removing it manually?
I tried using geometry checker


texts:
Is there a way to get the DE-9IM matrix between two geometries in a Spatialite function, analogous to the st_relate() function in PostGIS?

(I find it very helpful to spatially join using intersect() and then do a count(*) group by with the DE-9IM matrix to determine what my interesections are like, whether on the edge or not.  Then I often calculate the ratio of areas for geometries that don't line up exactly to tell me about slivers.)

For reference:

http://en.wikipedia.org/wiki/DE-9IM

http://postgis.net/docs/ST_Relate.html


texts: Keep name of Local Authority when converting KML to GeoJSON?

texts:
I am trying to convert this KML to GeoJSON.

All I want to do is draw the map with D3 and as I hover over each local authority, have the name in a popup.

I tried this online tool and it seems to work, however I am loosing the name of places. Also when I later convert my GeoJSON to TopoJSON using this tool, it isn't a reduction in size of the 80% I was expecting.


texts: Saving QGIS Project file (*.qgs) with PyQGIS?

texts:
I'm having trouble with QGIS (2.18.12) project files on Windows, if saved via python command. The project file is saved via:

project = QgsProject.instance()
project.write()

The file is saved, I can open it and it contains all the layers it should.

However if I first use the save command (described above) and then press the QGIS 'Save button' (the disk symbol) I get an message saying:

'Project file was changed' ... (see image)

This confuses me as I'm in fact overwriting the existing project file. If I reopen the file it will show me the newly saved status.

Besides this the actual message doesn't seem to change anything. Whether I press 'OK' or 'Cancel' the actual file content of the project file, doesn't change or revert (if I press 'Cancel' and close). When reloaded the qgs-file always shows the correct status (= the project stats when the python save-command was executed).

The date/time in the message seems to refer to the date I originally saved the project file via the QGIS 'Save Button'.

If I press the save button again the message won't appear again (until I use the python save-command). Also, if I save the project via python, then close the project, reload it and press the QGIS 'Save Button' the message won't appear.

How can I save a QGIS project file and use the QGIS 'save button' afterwards without this message showing up?


texts: No file is created when I choose 'Save style ...' in layer properties

texts:
First of all, I create a layer using the Add Delimited Text Layer plugin. The layer contains longitudinal and latitudinal data and a number that is assigned to each data point. When using a graduated color to style my data points, I define classes based on a specific range. I'd like to use the same classes for all of my layers. Therefore, I want to save my style using the 'Save style ...' button, in order to be able to load it into the other layers. However, when saving the .qml file, no file is created in the specified folder (e.g. desktop).

Another strategy may be to set the current style as default by using 'Save As Default'. When clicking 'Restore Default Style' in the other layers, a prompt shows up that says 'style not found in database'.

Thanks for helping me out!


texts: Converting polygon layer to *.txt files using ArcGIS Desktop?

texts:
Is there any way to convert polygons that do not have coordinates in the attribute table to a txt file with coordinates.

I will be using the txt file in a different software (not ArcGIS), so I need the txt file with polygon coordinates.


texts: Adjust Picture Marker Symbols with Python? (ArcGIS 10 SP1)

texts:
I am working toward a map cache, and therefore have about 12 versions of the same datasets to be shown at different scales.  Invariably, the symbol sizes I have set on highway shields (for example) need to be fussed with in order to make them look right.  In order to change the symbol size (scale to text size doesn't fit my application), you have to click four times from the label manager, and you have to do it over and over and over until all eight label classes are done.  And then you want it just a tad larger....

Can python (or anything else) be used to adjust picture marker symbol size across multiple label classes?
Thank you!


texts: Using cursor in ArcPy gives TypeError: cannot concatenate 'str' and 'Polygon' objects?

texts:
I have been trying to use arcpy.da.seachcursor to rip through a datasets rows and use each in an intersect analysis like so:

with arcpy.da.SearchCursor(in_feature, ['SHAPE@', 'NAME']) as cursor:
    for f in cursor:
        data = []
        print '...Working on {}'.format(f[1])
        for path, subdirs, files in os.walk(in_folder):
            for i,x in enumerate(files):
                if x.endswith('.shp'):
                   arcpy.Intersect_analysis(os.path.join(path, x)+';' + f[0], 'in_memory\\r{}'.format(i))
                   data.append('in_memory\\r{}'.format(i))
        arcpy.Merge_management(data, out_folder + '\\out_dataset_{}.shp'.format(name))
        arcpy.Delete_management('in_memory\\')

it yields the following error code message:

TypeError: cannot concatenate 'str' and 'Polygon' objects

Upon reviewing how to fix this, I have stumbled upon code like:

row[0].overlaps(row2[0])

that is used by da cursors (reference: https://arcpy.wordpress.com/2012/02/01/find-overlapping-features/). This has piqued my interest but I cannot find documentation on the implementation of this data structure. Where do I find it? Are there methods for doing other spatial processes like intersect?

Hopefully through this I find out the answer to my error code, but if anyone has a quick bit of knowledge outside of the data structure I am looking for info on, please drop it :)


texts: How to get centerline of each polygon (rectangles) in QGIS

texts:
How to get centerline of each polygon (rectangles) as in Figure 1 in QGIS to obtain centerlines as in Figure 2

Figure 1

Figure 2


texts: Getting EPSG from Feature Classes of Geodatabase?

texts: Where can I find the EPSG code of a Feature Class?


texts: How to cartographically represent water flow?

texts:
I am modelling water flow on steep slopes for an erosion study.  You can see from the screenshot the results from the ArcGIS-based flow analysis.  To produce these flow lines, I generated a depressionless DEM, calculated flow direction and finally flow accumulation.  The raster flow accumulation product was vectorized with ArcScan.

Although this method does a great job at showing water paths, it does not provide any visual measure of magnitude or accumulation.  How can I better represent the flow of water across steep surfaces so that magnitude and/or accumulation are represented?


texts: Converting polyline to raster with varying cell value along shapes?

texts:
I have a polyline file with start and end elevations. I'm trying to convert this to a raster which will have as values a gradient between the start and end elevations along the line.

Is there any way to do this?


texts: r.sun GRASS GIS .. calculate sunhours for a specific coordinate pair

texts: I am wondering if I can calculate the sunhours for a specific coordinate pair with r.sun on GRASS GIS. I want to give a coordinate pair as input and calculate how many sunhours there are on a specific day.


texts: Can I use relative path for image in SLD 1.1.0

texts:
Can I use relative path in the SLD 1.1.0? When validating it in geoserver, it said

Parsing failed for ExternalGraphic: java.lang.IllegalArgumentException: URI is not absolute

SLD:

<se:PointSymbolizer>
      <se:Graphic>
       <se:ExternalGraphic>
         <se:OnlineResource xlink:type="simple" xlink:href="images/mrtstation.png"/>
        <se:Format>image/png</se:Format>
       </se:ExternalGraphic>
       <se:Size>4</se:Size>
      </se:Graphic>
</se:PointSymbolizer>


texts: How to put image on top of OpenLayers 3 map by pixels

texts:
I need to put image which should be in middle top of openlayers 3 map. So, it should be defined with map canvas and pixels. Is that possible to do? This is my code:

 var arrow = new ol.layer.Image({
        source: new ol.source.ImageStatic({
            imageExtent: bounds,
            url: imageUrl,
            imageSize: size

        })
    });
    map.addLayer(arrow);


texts: Creating a bar chart in QGIS 3 from a list attribute


texts: I have been searching for a watershed layer for California or for the entire US that encompasses all upstream area from ocean outlets (e.g. the entire catchment draining into each ocean output). I have found the different HUC levels in the Watershed Boundary Dataset, however, these either group multiple outlets together or do not represent all upstream drainage area. There are a few different ways that I can imagine going about making this layer, however, I thought I would ask here first to make sure I am not reinventing the wheel.


texts: Combine duplicate geometries into single records with v.clean in QGIS 3.6

texts:
I have a points file that contains multiple duplicate geometries.  I am using the QGIS "Processing Modeller" to clean up data, and want to identify duplicate point geometries and combine so that I have a single entity.

I used the "Delete Duplicate Geometries" command, and had a significant reduction in the number of records.  I was skeptical at first so used the GRASS command v.clean with "rmdupl" command and received two outputs:

Cleaned
Errors

I noted that the "Cleaned" file still had duplicates in it.  Is there another step to remove the "errors" from the "Cleaned" file that I am missing.

As it turns out the "Delete Duplicate Geometries" appears to have given the correct answer, but I am still unclear using v.clean as it looks like a powerful tool.


texts: Unequally spaced points along a line script - QGIS Python error trying to add integer attribute

texts:
I am trying to plot points at specific positions along a line (which exists in the 'current layer' and is selected on the canvas) - these positions are irregular, so I am intending to loop through a list (chainage) and then add a new point layer.

I have created a python script by scrounging around various sources. You may recognise the code as mashed from yours.  To get it to nearly work in QGIS 3.4 I've had to adapt some of the existing code to work with the changes to the API since some of the code was written... Introduced errors within this are mine. However, the problem appears to be that it cannot add the integer attribute, even though I used the QGIS pycookbook for that bit:

    from qgis.core import (QgsFeature, QgsGeometry,
                       QgsVectorLayer, QgsProject,
                       QgsField)
from qgis.PyQt.QtCore import QVariant
from qgis.utils import iface
#list of positions along line (chainage) - to do - read a csv file!
chainage = [0,13,34,54]
# Create a new memory layer and add a distance attribute
vl = QgsVectorLayer("Point", "distance nodes", "memory")
pr = vl.dataProvider()
pr.addAttributes( [ QgsField("distance", QVariant.Int) ] )
layer = iface.activeLayer()
selection = layer.selectedFeatures()
for feature in selection:
    geom = feature.geometry()
    # Loop though all the chainage distances
    for dist in chainage:
        fet = QgsFeature()
        fet.setGeometry(geom.interpolate(dist))
        fet.setAttributes(dist)
        pr.addFeatures([fet])
        vl.updateExtents()

QgsProject.instance().addMapLayer(vl)

This nearly works and then falls over at line 19 an error related to the integer.

  TypeError: QgsFeature.setAttributes(): argument 1 has unexpected type
  'int'

If I comment out that line it works, but obviously I don't have the distance as an attiribute, just the geometry.

Maybe it's because I've cherry picked from different places / versions, because I'm not a python programmer.

I've included all my code because it might be useful to others... There seem to be a few similar questions out there.

I know my code is clumsy and not 'best practice' - once I get it working I can prettify it and probably rationalise it a bit.


texts: Setting up Python 2.7 to work in both QGIS and ArcGIS Desktop?

texts:
I have both ArcGIS 10.1 and QGIS 1.8 installed and both include Python 2.7 in slightly different setups and in different folders (c:\Python27\ArcGIS10.1\ and c:\Program Files (x86)\Quantum GIS Lisboa\apps\Python27).

It is only ArcGIS Desktop which include python.exe, while it is really QGIS I would like to start experimenting with for some scripting, but preferably just running my scripts from the command console.

How do I do this?

I'm thinking of settings environment variables etc.

I am new to both QGIS and Python.


texts:
I am having problems with r.viewshed.cva in GRASS GIS 7.0.3 and 7.0.4. I run the addon with the DEM and 100 random vector points as input. The vector points are composed by the cat, and X and Y colum with coordinates. Finally I get this error:

Traceback (most recent call last):
  File "C:\Users\AppData\Roaming\GRASS7\addons/script
s/r.viewshed.cva.py", line 174, in
    main()
  File "C:\Users\AppData\Roaming\GRASS7\addons/script
s/r.viewshed.cva.py", line 140, in main
    if flags['k'] and options["name_col"] is not '':
KeyError: 'name_col'

My laptop is a Core i7, Windows 10.


texts: Finding DEM for Greece (in high resolution) in order to calculate road slopes?

texts: I want to calculate the slope of road segments and I need dem for Greece. Where can I find it? I have found another one but the slopes were huge which is impossible.


texts: ArcCatalog database connection showing blank tables

texts:
I received info to connect to a database (without a spatial component) that has tables I need for a project I'm working on. Problem is once I establish the connection in ArcCatalog I can see all of the tables contained in the database, but when I try to open all the fields are there but they aren't populated with any information. The person who created them says it looks good on his end.

I attempted to export one of the tables to a .dbf file and it gave me a 999999 error with a 000592 warning.


texts:
i want to show my shapefile on the web
my template file is:

<title>Simple Mapserver Template HTML Document</title>

my map code is:

MAP
NAME "parcels_map"
STATUS ON
SIZE 400 300
EXTENT 489885 5450946 490904 5451637
UNITS METERS
SHAPEPATH "shpfile"
IMAGECOLOR 255 255 255

WEB
TEMPLATE "template.html"
IMAGEPATH "C:\\wamp64\\www\\parcel\\output\\"
IMAGEURL "/output/"
END

LAYER
NAME "parcel_boundaries"
DATA parcels
STATUS DEFAULT
TYPE LINE
CLASS
NAME "parcels_color"
COLOR 255 215 0
END
END

# end of map file
END

i send a request to mapserv.exe in CGI mode and then my output is a broken PNG in my browser while a png file create in my output folder.
what should i do to solve my problem
please help me.I'm a beginner to mapserver and totally have disappointed


texts:
I have two questions regarding picture marker symbol in flex viewer:

How can I make a composite symbol so that users can re-size and drag it. Any example code could help a lot.
I have three symbols (2 text and a picture marker symbol) in a composite symbol. A picture marker symbol is drawn where the user exactly clicks on the map. One text Symbol on the right and the problem is with the left text symbol. I am not able to place it on the left side.
By setting its x offset property to 0, it centers it self in the picture marker symbol.
By setting x offset to negative the result is same as setting it to 0.

Any ideas?


texts: Importing submodule with python for QGIS?

texts:
So I am probably missing something really obvious here, but for some reason, I can't import submodules in my python plugin.

This is the general structure of my plugin:

├── __init__.py
├── metadata.txt
├── pluginMain.py
├── plugin_ui.py
├── plugin_ui.ui
├── resources.qrc
├── resources_rc.py
├── sub
│   ├── __init__.py
│   └── qcf.py
├── ui.py

I'm trying to import the qcf file from the sub module, but for some reason, this simple snippet fails to execute when QGIS calls it:

import sub.qcf
print sub.qcf

This is (part of) the stacktrace I receive:

File "/Users/.../.qgis2/python/plugins/.../ui.py", line 26, in
    print sub.qcf
AttributeError: 'module' object has no attribute 'qcf'

As you can see, importing works fine, but actually using the submodule does not work.

When I try to run the same code from within the python shell (in the top level directory), everything works as it should:

Python 2.7.11 (default, Dec  7 2015, 12:36:04)
[GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.0.72)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import sub.qcf
>>> print sub.qcf
<module 'sub.qcf' from 'sub/qcf.pyc'>

So, why does this simple import not work from within QGIS and is there any way to solve this?


texts: Identifying whether the user has clicked on a circle or a polygon?

texts:
I have a user-created graphics layer that might have lines, polygons and circles on it. When the user clicks to edit one of the graphics, I would like to be able to identify whether they have clicked on a circle or a polygon but the geometry.type for both of these is "esriGeometryPolygon".

Is there a straightforward way to tell the difference?


texts: I just installed Qgis2threejs 2.2 for QGIS 3.2.3 and it works fine but I can't see the options of the plugin anywhere. There is supposed to be a blue one but nothing.


texts: Some features in advanced digitizing tool are not active

texts: I'm working on QGIS 3.2.1 and when I tried to edit a layer using advanced digitizing tool and set the layer in the editing mode, the split features is the only active icon and the rest is not active as in the below screen. How can I change this?


texts: Downloading large data from PostgreSQL DB using WFS request in Shape-zip format using GeoServer

texts:
I am using GeoServer 2.18.0 on CentOS 7 & it's deployed under Tomcat 9.

We have few large spatial tables in our PostgreSQL/PostGIS database (those table has spatial index and vacuumed)
Using WFS request end-users downloading data as 'Shape-zip format'.
If file size reaches up to 2 GB then request gets crashed tried to find its max. limit & based on [geosolution site][1] article max. limit is 2 GB.

My Questions :

How to handle large data requests in GeoServer for 'shape-zip' format ?
Can we increase max. limit for shapefile ?
If not then what's the alternate options ? (thinking to use GeoJSON format)

[1]: https://geoserver.geo-solutions.it/educational/en/vector_data/otherformats.html#:~:text=Only%20one%20geometry%20type%20per,should%20the%20limit%20be%20exceeded)


texts: Using own basemap with ArcGIS API for Javascript?

texts:
I'm trying to make a map using the Javascript API.

I want to use my own basemap, but I can't find the most simple code to do this.

I always see something like this:

  function init() {
    map = new esri.Map("mapDiv", {
      basemap: "satellite",
      center: [-97.395, 37.537],
      zoom: 11
    });

But I don't want to use the ArcGIS online basemaps.

What is the code to make a map service I am hosting my basemap?


texts: ArcMap Append (Data Management) Tool is not working correctly

texts:
Software version: ArcMap 10.5

I am attempting to append data from a Shapefile into an existing SDE Geodatabase feature class using the Append (Data Management) tool. The tool runs, and the results say that it has been successful and no errors are shown, however the point geometry and attributes are not shown in the feature class.

I have tried in both an edit session and not, and using NO_TEST to map the attributes. I cannot see any reason why it wouldn't work.

What am I doing wrong and are there any alternatives to this tool if it is a known bug?


texts:
I'm using ArcGIS 10.3.1, and I want to create a mosaic dataset, but my rasters are in 12 bits, and a mosaic doesn't have this option.

What can I do?


texts: i want to load 2 layers from geonetwork in one basemap on geonetwork !! so my question turn about how to configurate the file config-viewer.xml to support 2 arguments name  in the same ows-context:Layer


texts: Computing length of network inside rectangle in PyQGIS?

texts:
I have a QgsRectangle and a network layer of lines. How is it possible to calculate the length (or percentage) of the network that lies inside the rectangle?

I found this answer that helps computing the percentage of a line inside a polygon, but I would love to avoid having to convert each rectangle in a polygon (I have almost 50000 rectangles) and iterating over that and over the whole line network (which is also some 25000 features).

Here is what I have:

from qgis.PyQt.QtWidgets import *
from qgis.PyQt.QtCore import *
from qgis.core import *
from qgis.gui import *
import itertools

registry = QgsMapLayerRegistry.instance()
network = registry.mapLayersByName('network')[0]

#all of this is to create the rectangles that I'll need to compute later

ext = network.extent()

x1 = ext.xMinimum()
x2 = ext.xMaximum()
y1 = ext.yMinimum()
y2 = ext.yMaximum()

X, Y = range(x1,x2,100), range(y1,y2,100)
XY = itertools.product(X,Y)

rectDic = {}
for x,y in XY:
    p1 = QgsPoint(x,y)
    p2 = QgsPoint(x + 1000, y + 1000)
    r = QgsRectangle(p1,p2)
    rectDic[(x,y)] = r

#and this is where I get stuck

networkPerRect = {}
for pt in rectDic:
    r = recdic[pt]
    lenOfNetwork = length of network inside r #NEED HELP HERE
    networkPerRect[pt] = lenOfNetwork


texts:
I have imported a Vector map into GRASS GIS. There are some problems with the topology and some polygon areas overlap.

I used 'v.clean' to break overlapping polygons and now I have multiple polygons mapped to the same category in the attribute table.

Is there any way to convert this feature to category mapping so that each feature maps to a separate category? I don`t care that I will have duplicate information in the attribute table.

All I want, is that starting from the original vector map (that has 1000 features and 1000 corresponding records in the attribute table), that after running 'v.clean' and obtaining 1500 features to have 1500 corresponding records in the attribute table.

Can this be done?


texts:
I am interested in mapping surface water bodies in a desert. I have used Tasseled Cap Wetness (TCW) in GEE. After running my code, I found that water bodies (locally called tobas) were not marked. I guess my choice of data may not be correct, my visualisation parameters may be wrong, may be it did not set threshold, or I chose a bad index.

var tobas = ee.FeatureCollection("users/tariqabbas/tobasId"),
    roi = /* color: #98ff00 */ee.Geometry.Point([71.46225179238947, 28.783863923754982]);

var TCAP = function(image)
{
var wetness = image.expression(
'(blue * 0.1511) + (green * 0.1973) + (red * 0.3283) + (nir *0.3407) + (swir1 * -0.7117) + (swir2 * -0.4559)',
{
  blue : image.select('B2'),
  green : image.select('B3'),
  red : image.select('B4'),
  nir : image.select('B5'),
  swir1 : image.select('B6'),
  swir2 : image.select('B7')
});
image = image.addBands(wetness.rename('wet'));
return image ;
};

var OLI = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterDate('2019-06-01', '2019-06-30')
  .filterBounds(roi)
  .filterMetadata('CLOUD_COVER','less_than',10)
  .map(TCAP);

Map.addLayer(OLI.median(),{min:0,max:1,bands:('wet')});
Map.addLayer(roi);


texts:
I need to obtain ground truth data of Google Earth for a certain region for the purpose of constructing confusion matrix. The data should correspond to Landcover classes.

My main aim is to construct a confusion matrix with rows and columns classified as land covers. Now I already have a experimental data set of latitudes and longitudes corresponding to each land cover classes. So I suppose I have my rows are filled. Now to obtain column values I need to fill it with ground truth data of Google Earth. However I am not exactly sure what does it mean in this scenario. I am guessing that I need to enter the lat-long values in GE and obtain land cover classes from it corresponding to the imported lat-long values.

Am I right? If yes then how do I get the land cover classes of the imported lat-long regions?


texts: I'm looking for spatial clustering algorithm for using it within PostGIS-enabled database for point features. I'm going to write plpgsql function that takes distance between points within the same cluster as input. At the output function returns array of clusters. The most obvious solution is to build buffer zones specified distance around the feature and search for features into this buffer. If such features exist then continue to build a buffer around them, etc. If such features not exist that means cluster building is completed. Maybe there are some clever solutions?


texts: Is it possible to create a default workspace in QGIS 2.0?

texts:
I have been using QGIS for quite some time. It always happens that I change the position of toolbars and panels to make the workflow more efficient. But then when I close QGIS and I have to start a new project, the position of toolbars and panels remain the same. I have to put them back individually to their position.

My question is, is it possible to create a default workspace? I mean, if I want to put back at once all the tools and panels in their position, I just have to click one or two buttons. Like in photoshop, they have several workspaces.


texts: Setting up ArcPy virtual environment with Pandas?

texts:
Do I need to bounce back and forth between two virtual environments in order to run Pandas and ArcPy in the same script, while ensuring Pandas (or accidental updates to NumPy) doesn't break Arcpy?

I'm getting a new computer, and I want to set up my Python environments correctly from the start. As an example, below is a common script that I normally write/run in PyCharm. (Somehow I've managed to not break ArcPy by running everything in "C:\Python27\ArcGIS10.6" but I know that is wrong, so I am trying to fix it.)

import pandas as pd
import arcpy

myxls = r'C:\Users\name\path\myfile.xlsx'

tab = 'Sheet1'
xl = pd.ExcelFile(pts_in)
df = xl.parse(tab)
# do things here
out_csv = r'C:\Users\name\path\mynewfile.csv'
df.to_csv(out_csv, index=False, encoding='utf-8')

working = r'C:\Users\name\path\myshapefile.shp'
arcpy.MakeXYEventLayer_management(out_csv, "Longitude", "Latitude", "pts_layer",
                                  spatial_reference=arcpy.SpatialReference(4326))
# then do other things


texts: In ArcGIS 10.2 I created two buffers using the default settings one with a distance of 1 kilometer. When I calculated the area it has a square kilometer area of 3.14. All seems fine there to me. However, when I run the calculate area for the buffer that is 2 kilometers the area is 1.25581. Interestingly, as you go up in size to 3 kilometers, etc. the area calculations continue to be one decimal off, that is the area of a circle with a radius that is 2km should be 12.5581, 3km should be 28.27 rather than 1.256 and 2.287 respectively. The projections are the same-albers equal area. Any ideas what is happening??


texts: Removing annotation layer according to data source using ArcPy gives NameError?

texts:
I work with this FGB annotation (ArcView 10.3):

and i try to remove it with this code:

import arcpy,os
from arcpy import env

env.workspace = r"G:\PROJECTS\road611\gis\aa"
for mxdname in arcpy.ListFiles("*.mxd"):
   print mxdname
   fullPath = env.workspace + "\\" + mxdname
   if 'env_70-42' not in mxdname:
      mxd = arcpy.mapping.MapDocument(fullPath)
      df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
      for lyr in arcpy.mapping.ListLayers(mxd, "*",df):
         if lyr.isGroupLayer == True: continue
         if lyr.dataSource == r"G:\PROJECTS\road611\gis\lyr\18_8_16\New File Geodatabase.gdb\anno":
            arcpy.mapping.RemoveLayer(df, lyr)
            print 'RemoveLayer'
      mxd.save()
      counter = counter + 1
del mxd

but i failed to remove the layer and i get this error:

NameError: The attribute 'dataSource' is not supported on this instance of Layer

but if i remove manually the "anno" feature class, and try to remove other FGD feature class, then i succeed and get this:

>>>
antiquities_70-42.mxd
RemoveLayer
1


texts: Finding duplicated Parcel according to two fields?

texts:
I have parcel layer (geodatabase), which has total number of 250,000 parcels, and I want check if the parcel is  duplicated.

In other words I want to find out the the rows which have the same Block number, Parcel Number in a certain  Location.

I am wondering if there is a method to do that?


texts: If road lengths are raster layer, how to calculate road lengths for different buffers?

texts:
A question related calculation of road layer lengths using raster layer for land use regression.

If raster layers are used instead of vector layer for calculation of land use area and road lengths, how can we represent the road lengths layer as a raster layer and how can we calculate road lengths in raster layer for different buffers?

How can we build a model for this using ModelBuilder?


texts:
I have a vector layer in which I have added some fields of type 'float' using arcpy.AddField_management. I'm performing simple calculations between those fields using update cursors but the problem is that when I divide the values by 0, I get a run-time error even though I'm using Inf with the float function.

This is my code:

with arcpy.da.UpdateCursor('flurstcuke_cropped5', ['FID','totresgeb', 'totresgbpt','difresarea','scaledif']) as cursor6:

for row6 in cursor6:

   row6[3] = row6[2] - row6[1]

   try:

       row6[4] = ((row6[2] - row6[1])/row6[1]) * 100

   except ZeroDivisionError:

      row6[4] = float('INF')

   cursor6.updateRow([row6[0], row6[1], row6[2], row6[3], row6[4]])

The message I get is that 'the value type is incompatible with the field type. [scaledif].' And I know it's referring to Inf since I have also set up a counter in my script to determine where the program is. But I don't understand that since 'Inf' is a floating number here and the field's type is also float


texts: Tracking views of Story Maps on ArcGIS Online?

texts:
Is there a way to track how much traffic (how many views) an ArcGIS Online (by Esri) story map gets?

Google analytics type thing?


texts: When pressing field calculator icon Ubuntu log out!

texts:
I'm having a weird issue, as said in the title when I press the field calculator button the operative system log out.
Anyone as experienced this issue?

I'm running qgis 2.18 on ubuntu 16.04


texts:
Have a colorful wms-layer and a polygon. Want to create layer-group for a wms-service. Under the polygon wms-layer should be gray-scaled. Did some tests with blending (VendorOption) But all I got was negative effect or some weird color-changing. But no gray-scaling.
Alternative to make polygon transparent-gray looks awful, too.
Any suggestions?


texts: QGIS Least Cost Distance Points to Destinations via Roads

texts:
I have a layer of starting points (people's residences) and a layer of destination points (stores). I also have a roads layer with a cost field calculated from distance and speed of each segment.

What I'd like to end up with is a single layer, or multiple layers, with a cost field for each starting point to each destination point.

So if there are three starting points and three destination points, there should be 9 records in the output. Ultimately, there doesn't need to be any geometry in the output layer, it could just be a csv file. In reality I'm not even sure what type of geometry would be possible as an output.

This could also be run as a batch process, with each starting point as a separate layer. Just a thought if it simplifies things.

I've tried GRASS GIS v.net.distance but I had a bit of trouble with it. A few times it didn't execute properly, and once I ended up with a single line output layer, connecting the starting points to the closest destination. I think this tool does exactly that; finds the closest to feature from each from feature.


texts: Reprojecting vector layer from EPSG 3857 to 4326?

texts:
Based on examples from:

http://acanimal.github.io/thebookofopenlayers3/chapter03_10_reading_writing_features.html

http://acanimal.github.io/thebookofopenlayers3/chapter07_07_editing_features.html

I have managed to create simple editor, which pastes me created vectors to the textarea and then I created buttons, which downloads this GeoJSON string in a file.

I don't know why, but my app still gives me created vectors in EPSG:3857 instead of EPSG:4326 and it doesn't care that I defined EPSG to 4326 in both vectorLayer and for View.

var vectorLayer = new ol.layer.Vector({
            source: new ol.source.StaticVector({
                format: format,
                **projection: 'EPSG:3857'**
            })

´

var map = new ol.Map({
    ...
    })
  ]),
  view: new ol.View({
    **projection: 'EPSG:4326',**
  })
});

The problem is, that common user isn't clever enough to work with GeoJSON (download, open in QGIS, set up EPSG:3857 and save as EPSG:4326 in some other format), so I decided to create second button (for downloading) and defined var format = new ol.format.KML(); which is easy to use (double click opens and shows it in Google Earth), but this is not possible, because this works only with WGS84 coordinates in that KML file and not with other systems.

Is there (in OpenLayers) any possibility to easily create new duplicate layer with different coordinate system (geometry will be reprojected and I can just export)?

For example something like this could be nice:

var DuplicateWithDifferentEPSG = ol.proj.transformVectorlayer(vectorlayer, 'EPSG:3857', 'EPSG:4326')

How to simply reproject everything in my layer?

Any ideas?

Should I post exact copy of my code?

I don't get any errors.


texts:
I am trying to convert an imagine file with meter vertical and horizontal units to a foot projection.  Using gdalwarp, I am able to convert the horizontal but the vertical units stays in meters.  I did not see an option to either convert the z values or put the vertical unit type in the header.

I noticed gdal_grid has a -z_multiply option, but I would rather not have another step if possible.

Here is an example of my gdalwarp command:

gdalwarp -t_srs ESRI::C:\state_plane_foot.prj -tr 3.3 3.3 -r average -of HFA C:\utm_meter.img C:\state_plane_foot.img


texts: Why Landsat ETM+ panchromatic band' wavelength steps over visible range?

texts:
The Landsat ETM+ band-8 (Panchromatic) is identical to Landsat-8 OLI' Panchromatic band-8 in terms of spatial resolution i.e. 15 x 15m cell size. However, there is great difference in the wavelengths of both bands; ETM+ .52 - .90 and OLI 0.503 - 0.676 (micrometers).

See https://landsat.usgs.gov/what-are-band-designations-landsat-satellites

Clearly, the wavelength placed for ETM+ steps over visible range. A visual comparison between both bands also indicates outcome of these differences.

Obviously, it is noticed that OLI' pan band is much handy in visual interpretation as well as suitable for pan-sharpening and image classification.

There might be some good aspects behind ETM+ panchromatic wavelength extending beyond visible range, I am interested to be enlightened about the same reason.


texts: Spatial join of building footprints with Shapely and Fiona

texts:
I would like to do a spatial merge of the address nodes of Baton Rouge into their surrounding building footprints. Every building (polygon) should get the properties of the address (point).

Here is an visualization of the task:

I've tried to use the shapely+fiona solution from here. I've only changed it that the polygon shapefile get's written with the point shapefile's properties (before it was the other way round).

import fiona
from shapely.geometry import shape
from copy import deepcopy

with fiona.open("Building_Footprint.shp", "r") as buildings:

    with fiona.open("Street_Address.shp", "r") as addresses:

        outSchema = deepcopy(buildings.schema)
        outSchema['properties'].update(addresses.schema['properties'])
        print(outSchema)

        with fiona.open ("Joined.shp", "w", buildings.driver, outSchema, buildings.crs) as output:

            for address in addresses:
                for building in buildings:
                    # check if point is in polygon and set attribute
                    if shape(address['geometry']).within(shape(building['geometry'])):
                        print("yes")
                        building['properties'] = address['properties']
                    # write out
                        output.write({
                            'properties': building['properties'],
                            'geometry': building['geometry']
                        })

But I'm getting ValueError: Record does not match collection schema: odict_keys(['ID', ...]

Why is it failing?


texts: Open a KML with MultiGeometry's in QGIS 3.12

texts:
I tried to open the Sentinel-2 tiling grid (KML 100MB) from ESA's website using the Data Source Manager (Vector) in QGIS 3.12. There is no error, but also the geometries are not displayed. The file is opened and layers appear in the layers panel:

The KML has MultiGeometrys, which according to this answer cannot be opened in QGIS. How can I open the geometries found in this KML?


texts: How to use ST_SymDifference for polygons from the same table having the same name value in CARTO?

texts:
I am very new to PostGIS2.4 and Postgresql 10 and having hard times to make holes inside a polygon using another polygon from a single table liveuamap by using ST_SymDifference.

The polygon in the middle should be cleared.

I searched a lot for an answer but nothing could serve in my situation for CARTO. I used the following query :

SELECT a.cartodb_id, b.cartodb_id,
    CAST (ST_AsText (ST_SymDifference (a.the_geom, b.the_geom) )
    AS VARCHAR(350) ) SYM_DIFF
  FROM liveuamap a, liveuamap b
  WHERE a.cartodb_id = 105 AND b.cartodb_id = 116


texts: Detecting densely labelled areas in ArcPy?

texts:
I'm trying to detect densely labelled areas  for a given scale like explained below. so dense labelled area is the area having high rate of overlapping labels.

I'd like to detect densely labelled area to make inset maps to make labels clear for that areas. like explained in Map automation, create and populate inset maps for densely-labelled areas

I'm using ArcPy.


texts: How to split line into specific number of parts?

texts:
I've seen a lot of questions involving splitting a line with the help of a point layer.

I want to split a line into fractions of its length.

For example, I have a line 400 meters long, I want to split it into four lines of 100 meters long each.

There is the grass module v.split, but I get an error message when I start it from the qgis toolbox :

*"TypeError: object of type 'NoneType' has no len()"*

So I'm not sure if I get it to work if this would be a solution.


texts:
I get a random number of errors that come back as NA when I try to geocode a lot of places using rgooglemaps getGeoCode function. Can anyone tell me why? (Reproducible code below)

library(RgoogleMaps)
library(foreach)
###Replicating a large search data###
PlaceVector <- c(rep("Anchorage,Alaska", 20), rep("Baltimore,Maryland", 20),
rep("Birmingham,Alabama", 20))
iters <- length(PlaceVector)
###Looping to get each geocode###
geoadd <- foreach(a=1:iters, .combine=rbind) %do% {
  getGeoCode(paste(PlaceVector[a]))
}
geoadd <- as.data.frame(geoadd)
geoadd$Place <- PlaceVector

I get a random number of errors, usually around 15 where the latitude and longitudes in data frame geoadd come back as NA. I could loop it back on the NA's but that seems utterly inefficient. Do others have the same problem with the sample code provided?


texts: Semi-automatic classification plugin (SCP) error: local variable 'u' referenced before assignment

texts:
I use Ubuntu Xenial operating system, I can't use the Semi-automatic classification plugin (SCP) in QGIS. When I use a tool in the complement QGIS it just shuts down and then I need to restart the software. So I reinstall the complement then the software show me this message from Python:

local variable 'u' referenced before assignment

any idea of how can I solve it?


texts: ArcGIS Pro Network Analyst Routes appear to be incorrect?

texts: Using Network Analyst on ArcGIS Pro and a preconfigured streets network, I am trying to use the "Route" tool to solve the fastest way to service 20 different routes with 7-20 locations per route. However, each route I have solved goes completely against intuition. The stops are sequenced in ways that backtrack and cost more time than they should. For Sequence in the travel settings, I am using "Find Best." For input locations, I have separate feature classes for each of the twenty layers that contain all of the points in that route. I am trying to figure out if I am missing a crucial Network Analyst setting or if the network layer I am using might have been configured incorrectly. I included a picture of the shortest route just to showcase the type of problem I'm having.


texts:
I am currently following a tutorial book that is a bit outdated, therefore I have looked at the openlayers website:
https://openlayers.org/en/latest/examples/wms-tiled.html to find how to make a tiled WMS map.

I have had a problem with the imports that was already answered in another question on gis.stackexchange:
Unexpected token Identifier error - Import OpenLayers

I am looking to reproduce the behaviour on the website. However no matter what I do now I just get a blank screen. Even if I use the copy from website function and paste that into my own example.

I am afraid I still do not know enough about OpenLayers and mapping in general to fix this on my own.


texts: QGIS 3.16 Rasterize doesn't show anything

texts:
To give more info. The input layer is bus stop buffer dissolved which i buffered several stops in the district for 200 m and dissolved it to get those greenish polygons in the layer.
The output extent is the cip muara layer (in yellow, district shapefile layer). After running it i got the raster layer "rasterized" with... nothing shown on the layer. Did i get something wrong?


texts: Spacing in attribute table of ArcGIS Pro?

texts: For some reason when I import data from ArcMap to ArcGIS Pro, some of the data got messed up. For example; In ArcMap a column has "Concrete Foundation" ....but once imported into Pro, it ends up saying "ConcreteFoundation" and so on. I really need to learn Python apparently...but how can I fix this spacing error at the moment since I don't necessarily have the time to be studying Python at work.


texts: GDAL virtual file system, how to create your own

texts:
I just found a blob post about, GDAL virtual file system, and I was wondering if I could create my own and how to do so?

I want to create a virtual file system for reading directly from blob storage. So any example with virtual file system, i think I can managed to learn and create my own.

Do anyone have examples. (its microsoft azure I am going to create for, maybe someone already created one for aws s3)


texts: Import MGRS vector (old projection) to QGIS

texts:
I have large Excel tables of MGRS coordinates. These data were recorded from grid references based on maps  created in India 1960 projection. Is there any way to import these into QGIS to a) visualise without needing to convert to some X, Y / lat / long systems and b) reproject to WGS84?

Bear in mind that MGRS is self contained within just one field (ie, not the typical double X, Y fields)


texts: PostgreSQL is refusing connections from MapServer

texts:
I have MapServer 6.04 running in CentOS 6.7.  I am able to serve shapefiles and GeoJSON, but when I try to connect to a PostGIS (2.1.8)  enabled PostgreSQL (9.4.5) database it fails with the error: Database connection failed. I can connect to the database via QGIS and python from other machines. I can also run ogrinfo PG:"user=radar dbname=radar" ktlx -summary and get the following output:

Layer name: ktlx
Geometry: 3D Multi Polygon
Feature Count: 6
Extent: (-98.807300, 32.737200) - (-93.488900, 36.542300)
Layer SRS WKT:
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
    SPHEROID["GRS 1980",6378137,298.257222101,
        AUTHORITY["EPSG","7019"]],
    TOWGS84[0,0,0,0,0,0,0],
    AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
    AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
    AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]]
FID Column = index
Geometry Column = geom
radar: String (0.0)
val: Integer (3.0)

The connection info in my .map file is:

CONNECTIONTYPE postgis
CONNECTION "host=127.0.0.1 dbname=radar user=radar"

Every reference I have found in the documentation indicates that the the CONNECTIONTYPE and CONNECTION strings I am using are correct.  Is there a different connection string to use?


texts: Editing End Points in Polyline using ArcPy?

texts:
I am trying to make a polyline endpoint snapping tool in Python using ArcGIS 10. using a Search Cursor I read in the values to an array of line objects (which contain the feature id of the line, and two "EndPoint" objects - first and last - which contain the X and Y coord of each), then process some data interpolation of the points to give all points their correct values.

So far, that part works perfectly. However, when I try to go back and edit the shapefile, I can see the data that I am changing, but I can't seem to edit it.

Here is the code for updating the fields:

# --------------------------
#
# Update the file with the new EPs
#
# --------------------------

# Create update cursor
#
rows = arcpy.UpdateCursor(outputDirectory + "\\" + trails_fc + ".shp")

# Enter for loop for each feature/row
#
i = 0
for row in rows:
    # Create the geometry object
    #
    feat = row.getValue(shapefieldname)

    partnum = 0

    if feat.getPart(partnum)[0].X != allLines[i].startEP.x:
        arcpy.AddMessage("Change: " + str(feat.getPart(partnum)[0].X) + " to " + str(allLines[i].startEP.x) )
        feat.getPart(partnum)[0].X = allLines[i].startEP.x

    rows.updateRow(row)

    arcpy.AddMessage("Now is: " + str(feat.getPart(partnum)[0].X))

    i+=1

My readout consists of statements like this:

Change: -105.512166832 to -105.699533165
Now is: -105.512166832

For some reason, the rows are not updating. And for the life of me, I can't find a tutorial or instructions on how to edit a particular point in a polyline. I can only find how to edit a point as a field in a point shapefile.

Anyone have any ideas?


texts: Migrating PostGIS polygons into topologies?

texts:
With PostGIS, will migrating polygons into topologies using toTopoGeom() provide any benefit with regard to point-in-polygon search?

I have a lot of data, roughly 15 sets of data. And, I am wondering if there would be a benefit to migrating all of it to topology. Currently, 100% of my queries are point-in-polygon. However, most of these queries could benefit from some kind of subdividing and simplification.

Take for instance the TIGER States file. If 100% of your workload included finding out what state a point was in, what would be the benefit of migrating those polygons into topologies? What I'm thinking is I can avoid overlaps in the event of simplification, and that syntactically I'm concerned only with which state a user in in.

I'm just not sure I understand topology and its benefit to polygons, but it's appealing to me to know that by defining the edge any form of movement of that edge results in maintaining two faces.

This is demonstrated in section 13.23. and 13.24 of PostGIS in Action, 2nd Edition by Regina O. Obe and Leo S. Hsu. The pictures above are taken from the text.


texts: Attaching Excel sheet with different language font(hindi) to shapefile

texts:
When I try to join a excel sheet with Hindi font to a shp file, I get '?????' as a result in those cells which has Hindi font. (Before/Aft images attached)

I tried changing the locale language to Hindi but still no relief.

I suppose I am missing the right procedure. Can anyone help me with this..


texts: Training signature file for supervised classification in ArcGIS Desktop?

texts:
I want to do supervised classification in ArcGIS 9.3. I have a 321RGB composite of a pretty good quality which I want to use for supervised classification.

How can I create signature file which contains training sites, and how do I create training sites?

I do know what qualities are needed regarding the areas and the number of areas (number, all possible categories, spectrally separable etc.)


texts: Cannot find NATO UTM location in Sentinel-2

texts:
Regard coordinates 31.96212, -103.004715

UTM converters give it's UTM coordinates are 13/R/FR.

Example converter is here: http://www.rcn.montana.edu/resources/converter.aspx

But there are many of them and they give similar answers for these coordinates.

Simultaneously, in Sentinel-2 dataset here http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/#tiles/13/R/

I can't find FR subdirectory.

In google this location is here:

And finding the same place in Sentinel image browser I see, that tile is different

which stands for 13/S/FR i.e. the same UTM and square, but different band.

How is this possible?

UPDATE

KML with Sentinel-2 tiles also reports S tile in given location

UPDATE 2

According to this picture

taken from here, the FR square is located half in S UTM zone and half in R zone. Obviously, most automatic converters assign this square to R zone, while Sentinel-2 accounts it for S zone.

Is there any truth here?

UPDATE 3

The simple Python code, taken from here https://gis.stackexchange.com/a/224994/32207

bandVals = "CDEFGHJKLMNPQRSTUVWXX"

lon = 31.96212
lat = -103.004715

zone = int(lat + 186.0) / 6

if (lon >= 84.0):
    band = 'Y' if (lat < 0.0) else 'Z'
elif (lon <= -80.0):
    band = 'A' if (lat < 0.0) else 'B'
else:
    band = bandVals[int(lon + 80.0) / 8]

print '{:02d}{:s}'.format(zone,band)

also returns 13R.

Is this error in Sentinel-2 data or what?


texts: Building GEOS / GDAL halts on embedding manifest files

texts:
I'm trying to build GEOS and GDAL from source on a Windows machine.  I have Visual C++ for Python 9.0 installed with Lib / Include / PATH environment variables established (as far as I am aware). But I encounter a fatal error when trying to build GEOS and GDAL from source. I am using the following call to build GEOS in the GEOS root folder:

nmake /f makefile.vc MSVC_VER=1500

It compiles fine then exits on:

Microsoft (R) Incremental Linker Version 9.00.30729.01

Creating library geos_i.lib and object geos_i.exp
if exist geos.dll.manifest mt -manifest geos.dll.manifest -outputresource:geos.dll;2
NMAKE: fatal error U1077: 'if' return code '0xc000005'

After this error occurs, the source folder contains geos.dll as well as geos.dll.manifest. And, if I were to run the command nmake /f makefile.vc MSVC_VER=1500 again, it would pass on this instance but fail on a future one with the same error. By continuing to execute the nmake command I can eventually get through building geos.

I seem to recall that trying to compile each step by hand I got a more descriptive error message from mt saying that geos.dll could not be found. But when navigating to the folder, it is clearly there. And running the mt command by hand works just fine.

Any thoughts? I'm pretty new to building / compiling on Windows.


texts: Generating JPEG tiles using gdal2tiles.py?

texts:
I am using gdal2tiles.py to generate tiles and packing it into MBTiles using mb-util.
gdal2tiles by default generates PNG tiles. When I convert PNG tiles to JPEG large black area around tiles are left (due to lack of transparency in JPEG):

http://www.imagemagick.org/discourse-server/viewtopic.php?f=1&t=24048

Is there a way to directly generate JPEG tiles in gdal2tiles.py ?

I am using GDAL 1.9.2 on Ubuntu 12.04.


texts: Inserting image in attribute table and unable to save edits operation attempted on empty geometry

texts:
I have been trying to associate/insert images in an attribute table. I have created a field called pictures and data type is Raster. ( I have used this resource video to do the task)

Once I load the picture in the attribute table, every thing is ok. I can see the picture in the attribute table. But once I want to save the edits, I am getting this error: unable to save edits operation attempted on empty geometry.

ArcGIS Desktop 10.1 (Advanced), Windows 8


texts: Automatically assigning an attribute to point features within a polygon in ArcGIS Pro

texts:
I'm doing some keyboard-monkey level cleanup of a GIS with a vast number of point feature classes, some 50,000 points in total. In addition, I have a polygon feature class that divides the project area into distinct districts.

The point features all have an attribute (DISTRICT) that is intended to identify the district in which they lay, and this is supposed to correspond directly to the district name attribute (DISTRICT_NUM) in the district polygon feature class.

However, this DISTRICT attribute has many null or incorrect values, and will have to be redone in it's entirety. What tool or combination of tools could I use to automatically fill the DISTRICT attribute in the point feature classes based upon the DISTRICT_NUM attribute of the polygon those features lie entirely within?

Some limitations: I don't have anyway of adding new feature classes to the overall GIS, so processes that create new output feature classes will not work. Nor can I add additional fields to either feature class. That aside I have every license level and nearly every extension for ArcGIS Pro and ArcGIS Desktop.


texts: Get access to attributes of selected features stored in AGOL hosted feature layer from within Web AppBuilder

texts:
When you are passing one or more selected features into a geoprocessing widget via the ellipsis (...) in Web AppBuilder, what input parameter type do you use so that you can get access to the attributes of the selected feature(s)?

I have an AGOL hosted feature layer inside my Web Appbuilder (Developer 2.4) app that contains the attributes I'd like to work with.  I would like to assign an attribute to a variable so that I can perform some action on it from within the geoprocessing service.  I understand you can use the select widget to select one or more features and then pass those selected features as input to a geoprocessing service via the geoprocessing widget.


texts:
I have gtfs files from UK and I would like to know if there is a way to generate shapes.txt from these data.

I've found a tool for ArcMap, but credit is needed to use it.


texts: GeoNode 2.6 - Remove map and layers thumbnail?

texts:
When set a layer thumbnail and 500 error occurs. I am not sure what's going on. The quick fix may be just to remove the thumbnails from the layer and maps list page layout. Somehow this layout page, panels.html couldn't be located on the ubuntu server.

Any ideas?


texts: How to print to a specific scale with QGIS?

texts: How can I print a map with QGIS to a specific scale? E.g. I have a  map of Austria (projection: Albers equal area) and want a printout on paper at 1: 3000000?


texts: Filtering LiDAR data with las2las by scan angle range?

texts:
I am using las2las to filter .las data and I would like to filter the data between a scan angle of -15:+15; however the algorithm tells me the following: ?

  error bad lexical cast: source type value could not be interpreted.

I using the following as my input:

las2las -i /home/user/input.las -o /home/user/output.las --keep-scan-angle -15-15

How do I filter between a negative and positive range? From -15 to +15?


texts: Delineate watershed area for detention ponds

texts:
My study area has a series of small detention 'ponds' in the corner of fields. These are designed to capture runoff from the upslope contributing area (not from a discrete channel) during rainfall events. Each 'pond' has an embankment on one side (approx 1-1.5m high) to stop water from continuing downslope and detain it.  I would like to calculate the area draining to each pond using 1m resolution LiDAR data.
My first problem is that the ponds themselves are technically sinks. Am I correct in thinking that using the 'Fill Sinks' tool in ArcGIS with the Z limit set to ~1m should remove all/most other sinks but leave the detention ponds unfilled?
My second problem occurs when using the 'Watershed' tool in ArcGIS. I am using a raster representing the approximate area of each pond as the 'Pour Point', but this results in an output with three values, and a small catchment area with straight edges which is clearly incorrect. Should a pour point only have one cell, or is something else going wrong?
Many thanks.


texts:
Is there a correct way of including links to images within geoJSON?

This is the FeatureCollection I have for an object on my site:

{
  "type": "FeatureCollection",
  "features": [
    {
      "id": 1234,
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -0.322825,
          51.065361
        ]
      },
      "properties": {
        "popupContent": "In loving memory of\r\nGordon Wiseman\r\n15.4.30 ~ 11.8.06\r\nWho spent many happy hours\r\nin this park"
      }
    }
  ]
}

I also want to include the links to multiple images.  Is it OK to add another property? For example:

"properties": {
    "popupContent": "In loving memory of\r\nGordon Wiseman\r\n15.4.30 ~ 11.8.06\r\nWho spent many happy hours\r\nin this park",
    "media": [
        "https://example.com/image1.jpg",
        "https://example.com/image2.jpg"
    ]
}

Or is there a better / more standard way of doing it?


texts: Identifying south facing mountain land, above an elevation, within a slope range, less than fixed distance from major roadways

texts:
This odd question comes from believing that such an analysis is plausible in a GIS system.

Only, I don't how to accomplish the various discrete steps implied from the Title, nor where to get a topo map that already has the iso-elevation contour lines developed in some GIS format, nor how to color in areas which face SW-to-SE, nor how to decipher the physical slope from the proximity of iso-elevation lines.

All these things, sensible to the engineer in me, but wondering what set of stackexchange questions (or other resources) might even begin to answer my question....

...in my ideal world, I scan in a USGS Topo map for the region of my interest, and out the other end of my QGIS black-box procedure comes a map with polygons lit up as my areas of interest.

What tools and techniques....?  Learning resources....?

[and, if anyone wants to suggest 4 other appropriate tags, let me know....]


texts: How to select trenches only upstream of a station (select by location)

texts:
I have a shape (line) that represents the hydromorphology of watercourses  and contains several parameters. This morphological study of watercourses has been done on every 100 meters of surface watercourses.

In other side, I have another shape (point) with measurement stations that is superimposed on the watercourses.

For a specific analysis, I need 10 trenches (1000 m) in the upstream of each station to calculate the average of certain parameters.

I did a select by location with the method (are whithin a distance of the source layer feature), but it gives me a selection of 10 trenches in upstream and downstream of the measurement stations.

Anyone have a solution for having 10 trenches in upstream of each station only of the measurement stations?

The id's of trenches are not coded in a logical sequence.

The lines are connected to each other and flow in the source direction


texts: Leaflet see only one symbology at the time

texts:
I am working with Leaflet with nine layers and each have their own symbology like you can see in the image. For example, the first layer have 4 symbologys and you can see the colors in the map if you activate the layers, but you can see the colors in the legend all the time. I want that if you activate the layer, you can see the colors of the legend, but if i don't using this layer i only want to see the name of the layer.

Every layers is provided by a geoJson layer, i only find a solution for wms layer which is this   but there the symbology is an image and the layers are wms.


texts: Completely dicing feature class in ArcGIS Desktop?

texts:
Every so often I have to update my city's FEMA data, which requires running the Dice (Data Management) geoprocessing tool on the buffered floodway features so that it can be more easily processed as a feature class.

I run into the following problem every time I do it which requires a work around: when dicing the entire feature class, the output drops large swathes of features seemingly at random:

Input:

Diced output:

I am using a vertex limit of 500 when running Dice. Ive tried changing the vertex limit and this does not fix the problem, so i am forced to continue selecting the dropped portions and re running dice on them until the entire feature class is processed.

Anyone have experience with this issue or a possible alternate method?

I am running arcmap 10.3.1


texts: ArcGIS Spatial Join - "contains part of" match option

texts:
I would like to perform a spatial join using a (non-existent) "contains part of" match option. i.e. I want the join to include any and all polygons that overlap, not just touch at the edge. None of the match options available produces satisfactory results.

To clarify, in the example below, the yellow polygon should have a join count of 4 (A, B, C, D).

INTERSECT gives a join count of 6 (A - F)

The result I am interested in is the the join count per polygon

Any ideas for a work around or is there another tool that will do this?


texts: GeoServer - Single Image Request - Graphic Label Issues

texts:
I posted this on the OSGEO GeoServer Mailing List but received no love, so I'll try here.

I have GS 2.6.0. connected to a PostgreSQL 9.2 DB (PostGIS 2.1). I have an SLD where I am trying use graphic label options to label a bunch of weather conditions throughout 700 cities. I like how GS does not allow labels to conflict, but things are not working well. When I first open my OpenLayers 3 Map, the graphic label (icons) show nicely:

All the tiny black dots are cities I am try to overlay weather conditions. As I zoom in, I would expect different areas to be labeled, and when I zoom in ONE level, the graphics are still there.

But, when I zoom in one more level, and all subsequent levels, all the graphic labels disappear:

I am using OpenLayers 3.0 and am requesting a single image like this:

var openWeatherMap = new ol.layer.Image({
    visible: false,
    refresh: true,
    extent: mapExtent,
    source: new ol.source.ImageWMS({
            'crossOrigin':'anonymous',
            url: 'https://myserver.com/geoserver-2.6.0/myworkspace/wms',
            params: {'LAYERS': 'feed_open_weather_map'},
            serverType: 'geoserver'
    })

});

I've tried a million combinations of OL3 and SLD options. Some work better than other, but I can't ever get the graphic labels to appear beyond the zoom I mentioned.
My SLD looks like this: SLD Sorry for the pastebin. Does anyone know why GeoServer is not labeling with graphic properly?

Update: added getMap request samples.
The getMap request that shows the graphic labels is:

https://myserver.com/geoserver-2.6.0/myworkspace/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=feed_open_weather_map&CRS=EPSG%3A3857&STYLES=&WIDTH=1536&HEIGHT=1005&BBOX=-30056262.514183864%2C-19108413.37993557%2C30056262.514183864%2C20223023.89448472

When I put this in the browser, it shows the graphic labels.

The getMap request that does not show the graphic labels is:

https://myserver.com/geoserver-2.6.0/myworkspace/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=feed_open_weather_map&CRS=EPSG%3A3857&STYLES=&WIDTH=2880&HEIGHT=1005&BBOX=-12019026.271586256%2C102352.06492121238%2C16158719.83546112%2C9935211.383526284

When I put this in the browser, I only get the features (tiny dots) but no graphic labels:


texts: What is the best method to serve large GeoJSON data in OpenLayers?

texts: I would like to know about the best method to serve large GeoJSON data in OpenLayers. I am using extent-based loading and clearing the source when it is in out of extent. But still it is hanging while serving with big GeoJSON. I am accessing GeoJSON from GeoServer and ArcGIS Server based on extent.


texts: I downloaded about 9 LiDAR scenes from USGS and I want to use ArcGIS to merge them so that I can carry out tasks on them together as opposed to one by one.


texts: PyQGIS doesn't load user libraries properly

texts:
I'm trying to put my PyQGIS functions in separate .py files and load them from my main script.
from qgis.core import QgsProject

def getVectorLayerGeo():

    project = QgsProject.instance()
    layers = project.mapLayersByName('Yalcowinna')
    layer = layers[0]
    features=layer.getFeatures()

    coords = list()

    for feature in features:
        geom = feature.geometry()
        geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
        if geom.type() == QgsWkbTypes.PointGeometry:
            if geomSingleType:
                x = geom.asPoint()
                print("Point: ", x)
            else:
                x = geom.asMultiPoint()
                print("MultiPoint: ", x)
        elif geom.type() == QgsWkbTypes.PolygonGeometry:
            if geomSingleType:
                x = geom.asPolygon()
                print("Polygon: ", x, "Area: ", geom.area())
            else:
                x = geom.asMultiPolygon()
                print("MultiPolygon: ", x, "Area: ", geom.area())
        else:
            print("Unknown or invalid geometry")

    for coord in x[0][0]:
        coords.append([coord.x(),coord.y()])

    return coords

print(getVectorLayerGeo())

does exactly what it's supposed to, which is print a list of coordinates (as lists).
Then I try to call getVectorLayerGeo from another script:
import ee
from ee_plugin import Map
from qgis.core import QgsProject
from get_layer_coords import getVectorLayerGeo

vizParams = {'bands' : ['B4', 'B3', 'B2'], 'min':0, 'max':10000}
geom = ee.Geometry.Polygon(getVectorLayerGeom())
S2 = ee.ImageCollection("COPERNICUS/S2_SR")
image = S2.filterBounds(geo).first()

Map.addLayer(image, vizParams, 'Sentinel-2')
Map.centerObject(image, 10)

and I get
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 4, in <module>
  File "C:/PROGRA~1/QGIS 3.16/apps/qgis/./python\qgis\utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/Glen Shennan/Google Drive/African Lovegrass/ShannFlat\get_layer_coords.py", line 56, in <module>
    print(getVectorLayerGeo())
  File "C:/Users/Glen Shennan/Google Drive/African Lovegrass/ShannFlat\get_layer_coords.py", line 18, in getVectorLayerGeo
    geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
NameError: name 'QgsWkbTypes' is not defined

The thing is, when I copy getVectorLayerGeo into my main script it works as intended. Obviously I have the solution for this case but I'm intending to write a lot of functions and I don't want to put them all in one giant script.
QgsWkbTypes is part of qgis.core, but when I import qgis.core at the start of get_layer_coords.py script it doesn't seem to help. Am I doing something wrong or does QGIS not like calling functions from other .py files?


texts: IGlobeGraphicsLayer not updating in separate thread?

texts:
I'm currently using ArcGlobe to create a Polygon on a globe.  Using examples I've found through arcgis.com  I can now allow the user to click on the globe and create a polygon on the map.  This works as intended and the visualization on the map is seen.

However, I'm now trying to do this in a new thread so I can create a polygon overlay and attempt to return a list of Coordinates from a single point of entry (Since there doesn't seem to be a TrackPolygon function like there is in ArcMap).  When putting this code into it's own thread, the event handlers on the globe still work; however the GlobeGraphicsLayer.AddElement is no longer displaying this data on the globe.

Any ideas/suggestions?


texts:
I'm using QGIS 1.8 on windows XP machine and I try to do a concatenation of several fields to label the names of my local roads in the attribute table, but I can not find the correct syntax. The first field consists of om the street, the second type of link, eg: Avenue and the third orientation, for example, East or West. Could you tell me how to write the correct syntax using the dialog box of labels based on a formula.
  Attached a screenshot that shows the attribute table.


texts: Leaflet: obtain X and Y (pixel points) from GeoJSON lat, lon

texts:
I have a Geojson with Cities location. I'm trying to use the City Lat/Lon to send a WMS GetFeatureInfo to my Geoserver to add additional info to the label (temperature data), but I can't get the X and Y to pass in the request.
I've tried using the following code, but the ressponse is always different from a zoom level to another.
 var latlng = [GeoJsonLat,GeoJsonLon];

 var pixelPosition = map.latLngToLayerPoint(latlng);

Zoom Level 3:

Zoom Level 4:

And so on...
What am I doing wrong? Is there any alternative that I can use?
Update:
From the above code (lat, lon from GeoJson), pointing on Rome at ZoomLevel = 3, X=231and Y=340, while clicking on Rome with the mouse I get x=232 and y= 345

My Layer code:
cityLayer = L.geoJSON(data, {
  pointToLayer: (feature, latlng) => {
    let emptyGroup = L.layerGroup();

    let marker = () => {
      let tempValue;
      let html;
      let lat = feature.properties.latitude;
      let lon = feature.properties.longitude;
      var latLng = latlng;
      var point = map.latLngToLayerPoint(latlng);
      let bbox = map.getBounds().toBBoxString();
      var BBOX = map.getBounds().toBBoxString(); //coordinates
      var WIDTH = map.getSize().x;
      var HEIGHT = map.getSize().y;
      var X = point.x;
      var Y = point.y;
      console.log("X = ", X, " Y = ", Y);
      fetch(
        `https://www.xtreme-weather.com/maps/gfs/wms?` +
          "service=wms&version=1.1.0&request=GetFeatureInfo&layers=gfs:tempGfs&query_layers=gfs:tempGfs&styles=" +
          "&bbox=" + BBOX +
          "&feature_count=1" +
          "&height=" + HEIGHT +
          "&width=" + WIDTH +
          "&format=image%2Fpng&info_format=application%2Fjson&srs=EPSG:4326" +
          "&x=" + X +
          "&y=" + Y +
          "&time=" + validityTime +
          "&elevation=1002",
        { cache: "no-cache" }
      )
        .then((response) => response.json())
        .then((data) => {
          data.features.forEach((proerties) => {
            dataValue = proerties.properties.GRAY_INDEX;
            tempValue = dataValue.toFixed(2);
            // console.log(`<span class="my-div-span" style="color:white; width:100%;font-size:8px">${decode_utf8(feature.properties.name)}<br>${tempValue}</span>`)
            // return dataValue.toFixed(2)
            new L.Marker(latLng, {
              icon: new L.DivIcon({
                className: "div-icon",
                html: `${decode_utf8(
                  feature.properties.name
                )}
${tempValue}`,
              }),
              zIndexOffset: 10000,
            }).addTo(emptyGroup);
          });
        });
      // return tempValue
    };
    marker();
    // marker.addTo(emptyGroup)
    return emptyGroup;
  },
});


texts: How to add content below the map with Javascript?

texts:
Is it possible to enable some container below the map with JVectorMap? I have some links for all states which needs to appear when user clicks on state. I am new with Javascript so please help me.

For example, I have three links which needs to show when user clicks on Florida:
www.google.com
http://gis.stackexchange.com
ebay.com

How to define those links to be related to Florida? Links need to appear in the container when user clicks on that state.


texts: Set a flag as a user defined value on ArcGIS Portal

texts:
I have a geometric network published on a ArcGIS server where I'm trying to set a flag as a user defined value to perform a trace downstream task but I'm getting an error:

ERROR 000840: The value is not a Layer. Failed.

However if I'm setting the flag as a point layer the tracing works flawlessly.
Does anyone knows how can I tackle this problem?


texts:
I am using ArcGIS regression tools OLS and GWR for evaluating the association between dengue incidences and weather parameters, land use/cover, and demographic characteristics of the study area.

Poisson's regression is considered a better model for count variables as dengue cases. Can I perform the suggested analysis in ArcGIS?


texts: Changing legend in ArcGIS Online for use in a Dashboard

texts: Working in ArcGIS Online, I have a data set that has 3 types of line data in it , Plow, Grader, and Sweeping.  I have created copies of that layer and filtered out each type of the line work I want to show..... but in the legend it is still showing the other two sets of data.  How do I only show in the legend what I would like it to show...


texts:
I've followed this tutorial to create a custom Google Map. I've included a few other elements such as linking it up to Wordpress and clustering the markers.

It's all working great apart from the info in the info windows on each marker. I just can't seem to change the info within each one. I thought by changing the following lines it would change it but nothing affects it:

var html = "" + name + "
" + address;

This is the working map

Where can I put in my own custom data into the window? Also, if I could style the window on that would be even better.

It seems the clusterer is the problem, mainly this section, how can I take the html content and place it into the info window?

  // Change this depending on the name of your PHP file
  downloadUrl("<?php bloginfo('stylesheet_directory'); ?>/phpsqlajax_genxml.php ", function(data) {
    var xml = data.responseXML;
    var markers = xml.documentElement.getElementsByTagName("marker");
    for (var i = 0; i < markers.length; i++) {
      var name = markers[i].getAttribute("name");
      var address = markers[i].getAttribute("address");
      var type = markers[i].getAttribute("type");

      var offsetLat = markers[i].getAttribute("lat") * (Math.random() * (max - min) + min);
      var offsetLng = markers[i].getAttribute("lng") * (Math.random() * (max - min) + min);

      var point = new google.maps.LatLng(offsetLat, offsetLng);
      var html = "" + name + "
" + address;
      var icon = customIcons[type] || {};
      var marker = new google.maps.Marker({
        map: map,
        position: point,
        icon: icon.icon,
        shadow: icon.shadow
      });
      google.maps.event.addListener(marker, 'click', (function(marker, i) {
                    return function() {
                        infowindow.setContent(markers[i].getAttribute("name"));
                        infowindow.open(map, marker, html);
                    }
                })(marker, i));
      cluster.push(marker);
    }
    var mc = new MarkerClusterer(map,cluster);
  });
}

function bindInfoWindow(marker, map, infoWindow, html) {
  google.maps.event.addListener(marker, 'click', function() {
    infoWindow.setContent(html);
    infoWindow.open(map, marker);
  });
}

Specifically this, it's not putting the html content through the clusterer... at least this is actually changing the data in the window, just need to output the html content without breaking the clusterer:

      google.maps.event.addListener(marker, 'click', (function(marker, i) {
                    return function() {
                        infowindow.setContent(markers[i].getAttribute("name"));
                        infowindow.open(map, marker, html);
                    }
                })(marker, i));
      cluster.push(marker);


texts: How do I get the default NLCD symbology back after clipping/projecting?

texts:
NLCD landcover rasters downloaded from MRLC have default legends and symbology as shown here: http://www.mrlc.gov/nlcd11_leg.php

However, when I extract the raster through a mask, or reproject the raster, or do similar geoprocessing, the default symbology is replaced by random colors.

How do I get this symbology back? LANDFIRE offers a table with RGB values I can join to the file that, once joined, changes the raster image back to default colors when you go through properties -> symbology -> unique values -> add all values. Any similar solution for NLCD?


texts: Single numbers for accuracy and resolution for vector and raster data?

texts:
This question is meant as a modelling question, not having to do with particular data formats. It was edited to better separate spatial accuracy and resolution (precision).
In order to make the spatial accuracy and precision of spatial data in different datasets comparable, it would be very useful to be able to express spatial accuracy and spatial resolution as single numbers. That number ideally would be usable for both raster and vector data, so that data sets using both types of data become comparable by spatial resolution and/or by spatial accuracy.
When accuracy in each dimension is the same, the numbers for each dimension could be coalesced into a single number. But in some cases accuracy is significantly different between dimensions. A common example is vertical accuracy being different from horizontal accuracy. So perhaps for accuracy it is necessary to have the possibility of using a single number for each dimension.
Ideally, the unit of the numbers should be clear and also allow easy comparison of different datasets. For example, the unit could be specified in the coordinate reference system, or be an SI unit (to allow direct comparisons without needing to look up and convert data).
Can such definitions be found for both types of data that make the numbers comparable?


texts: Orfeo Toolbox for QGIS not activating after trying several suggestions

texts:
There are several other threads on this topic and after trying the suggestions I am still stuck with Orfeo Toolbox not activating in QGIS.
Running 2.16 right now... any chance updating to the newest 2.16.3 version would help?
Screenshot of my Orfeo's activate .


texts: What is the CRS for raster::area in R, and why divide by 1 million?

texts:
To calculate the area of polygon in R I have been using the following code (this is the example code from the package's documentation, mine is modified to contain my data):
p <- shapefile(system.file("external/lux.shp", package="raster"))
p$area <- round(area(p) / 10000000,1)
which returns
p$area
312.3 218.7 259.5  76.2 263.2 188.3 129.0 210.4 185.6 251.3 237.1 233.3
Anyway, these make sense as square km areas of the states of Luxembourg, and when I run this code with my data I get results that makes sense as well.
My question is: what is the CRS?
using crs(p) gives me CRS arguments: +proj=longlat +datum=WGS84 +no_defs  Warning message: In proj4string(x) : CRS object has comment, which is lost in output
From other posts I have learned that the warning code is nothing and that raster::area() returns square meters from lat long coordinates (such is the case here). But I was under the impression that WGS84 was a datum, not a CRS, so what is the CRS? Also why do I divide by 1 million? Is raster::area() converting to square millimeters?


texts: Difference spatial interpolation and trend surface analysis

texts:
What is the difference between spatial interpolation and trend surface analysis?

For example, I read a scientific study that applied a neural network to estimate a surface solely based on the x and y coordinates, claiming that they do spatial inteprolation.
However, to my understanding this is trend surface analysis, because Tobler's law is not considered.

Am I right?! Is the difference between these two that spatial interpolation takes neighboring observations into account?


texts: How to filter points based on proximity while inserting into PostGIS?

texts:
I have a large amount of locations that need to be inserted to a PostGIS database. However, I'd like to aggregate these points so that if there's already a location, say nearer than a kilometer, then to discard that point. I'd like a pointer or two on the appropriate query (or, if it can't be done in a reasonable way, then flame me...:)

I reckon it has to be done a bit like:

INSERT INTO locations VALUES (ST_GeographyFromText('POINT(54.55 26.33)'))
IF ST_Distance(ST_GeographyFromText('POINT(54.55 26.33)'), <all other points in "locations" should be evaluated here>) > 1000;

As most of the audience has probably understood by now, I'm not that bright when it comes to databases so don't be too harsh on me...


texts:
Does anyone know if you can limit what's shown in the GetCapabilities  document for a GeoServer WFS service?

I have restricted the workspace to read only mode (..r=* and ..w=NO_ONE). I hoped that would lead to that, for example CreateStoredQuery wouldn't show up in the GetCapabilities document, but it doesn't. I also only want to support GET requests, no POST requests.


texts: How can I setup a simple web mapping system?

texts:
I'm new to mapping and GIS and I'd like to setup a mapping system to display information on a map, I think the following could work but I'm not sure how good it is..

I store my data (lat, lng etc) on a PostgreSQL server
I build a GeoJson file containing the information I want to display on the map (using the data on the postgresql server)
I send the Geojson file to TileMill / PolyMaps / etc

Do I need to use tileStream or something to serve the tile ? To build on almost real-time system all I need is to rebuild the geoJson file right ?
is Wax useful for me ?


texts: Import and go through a large batch of Landsat images in QGIS?

texts:
I'm trying to reconstruct a glaciers frontal evolution through time.
To do so I use a whole lot (70+) of Landsat scenes which I got from USGS Earth Explorer.
I would like to visualize them as false colour composites and also with the NDSI (normalized-difference snow index).

What would be the most efficient way to get all these files from in the folders in QGIS and go through them through time?

EDIT: Currently, I load in all the bands that I want from a landsat scene and I do this for every single scene which is very time consuming. Then I give each band the colour that I would like it to have. Afterwards, I go through each scene one by one by checking on and off the bands that I want to see and draw a polyline of the glacier margin.


texts: Adding all Binary Rasters in a Directory Using Raster Calculator in a Standalone Python Script

texts:
I have a directory with binary raster files (1 band only) which I need to add in a standalone python script. The problem is the rasters have different number of rows and columns and different resolution (reading them in as numpy arrays and adding them does not work). When I use the raster calculator, it manages to do what I want, but now I need to use it outside QGIS and do it automatically for all rasters in the folder. How can I add all my rasters in my directory inputFolder and create one output raster? I have the following code so far:

inputFolder = "/home/usr/Desktop/"

def findRasters(path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join(root, file)

entries = []

for l in findRasters(inputFolder, '*.tif'):
    fileInfo = QFileInfo(l)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(l,baseName)
    QgsMapLayerRegistry.instance().addMapLayer(rlayer)
    layer = QgsRasterCalculatorEntry()
    layer.ref = l.name() + '@1'
    layer.raster = l
    layer.bandNumber = 1
    entries.append(layer)

expression = '(' + entries[1].ref + ' + ' + entries[0].ref + ')'
calc = QgsRasterCalculator(expression,'/home/usr/Desktop/final.tif',
                           'GTiff',
                           layers[0].extent(),
                           layers[0].width(),
                           layers[0].height(),
                           entries)

calc.processCalculation()

I am not sure how I should formulate my expression such that it takes into account all rasters in the folder


texts: Calculating the difference between two grids in percent with QGIS

texts:
I'm using the QGIS GDAL Raster Calculator to get the difference in percent between two grids.

Grid A contains an integer, and grid B contains an integer which may be the same as or larger than the corresponding cell in grid A. They are geographically identical (in terms of cell size, extent etc).

So in the raster calculator I use this formula:

((B-A)/A)*100

...expecting a percent value in the created raster showing the difference in percent between raster A and B.

But the created raster is almost completely zeroes, except for a few places containing 100% values (they do not correspond with my manual control samples). I've chosen Int32 as the datatype for the created raster (I don't need more precision than whole percent).

Is the formula incorrect for what I'm trying to do, or have I missed something else?


texts: Seeking QGIS version, or plug-in, for distributing QGZ project with simple user-interface and read-only user privileges

texts:
I would like to distribute read-only 3D mapping products to the non-GIS general public but am not sure of an easy way to do that with QGIS.  I see how to use Settings to create a customized user profile with limited functionality, but am not sure how to export/transfer this user profile along with my QGZ project file.  Many years ago I used an ArcView publisher software to create read-only GIS products consisting of pmf files and hidden database files.
Does QGIS have similar capabilities?
I'm currently using QGIS 3.10.6-A Coruna.


texts: Get lat, long from QGIS polyline shapefile

texts: I develop in C#. I want to develop a tool which needs to read the poly line coordinates. QGis is a very good tool that my users can use to generate a polyline .shp file. But how will I be able to access the coordinates(lat,long).


texts: Getting Projected bounds information in pyproj

texts:
I'm building an API where user can insert a number (3857, 27700,etc.) and my API will give the projected bounds.
I'm using pyproj's
crs_number = 27700
crs_custom = CRS.from_user_input(crs_number)
crs_custom.area_of_use.bounds

this function is giving me bounds but it is showing bounds like

-8.82,49.79,1.92,60.94

instead of

-90619.29,10097.13,612435.55,1234954.16

How do I get the expected bounds?


texts:
I was wondering if there are some free projects which enable extracting the building height information from openstreetmap data?

I downloaded data for a couple of large cities (Manhattan island, Munich...) and none of them contained the height information.

I googled a bit and found a project where Lidar point clouds have been overlayed on top of openstreetmap data, and then the height of the building rooftops have been identified by checking the point's Z coordinates.

Are there any other free projects which enable defining the height of the buildings of openstreetmap data? At least for world's major (largest) cities?


texts:
I have a geospatial python app that I want to use on Heroku.

This is my requirements.txt file:

Flask==0.10.1
Jinja2==2.7.2
MarkupSafe==0.19
Werkzeug==0.9.4
gunicorn==18.0
itsdangerous==0.23
numpy==1.8.0
requests==2.2.1
wsgiref==0.1.2

I want to use the heroku-buildpack-geodjango to install Python GDAL.
I create and do a git push.

heroku create --stack cedar --buildpack http://github.com/dulaccc/heroku-buildpack-geodjango/

git push heroku master

This all works fine. However the app crashes because it doesn't find GDAL.

These are my logs:

$ heroku logs
...
2014-03-23T18:57:08.224720+00:00 app[web.1]: ImportError: No module named osgeo
2014-03-23T18:57:08.348803+00:00 app[web.1]: 2014-03-23 18:57:08 [2] [INFO] Shutting down: Master
2014-03-23T18:57:08.348803+00:00 app[web.1]: 2014-03-23 18:57:08 [2] [INFO] Reason: Worker failed to boot.
...

How do I set up GDAL in order to work properly?


texts:
I’m new to working on automating processes with GIS data and I’m trying to come up with a solution in ArcMap that will automate merging polygons that have “Shape Area < 25” into larger polygons. I would like to merge it with the largest polygon that it is touching. I know dissolve is one solution to use to automate merging but in my case it would not work because the polygons that I’m working with are kind of in a grid pattern which needs to be retained. I would like to use a merge function similar to manually merging data via editor.

Does anyone know of a way using ModelBuilder?

Picture is an example so you can see what I’m trying to accomplish:


texts: How to chose which line goes on top of the other line

texts:
I am trying to do something very simple, but apparently impossible. Overlay two lines with different size and color signing same route different direction.
Like this one:
http://www.indianatrails.org/Pleasant_Ride/index.html

I couldn't find in CARTODB a way to put on top the green route.
Please help!


texts: Adding new fields to shapefile alters appearance of borders in QGIS?

texts:
I am using QGIS 2.6.0-Brighton and I am trying to add a custom field to my layer made from a shapefile of Italian provinces; everything seems fine when I create a categorized views using region codes, something like this:

then when I try to add a custom field to the associated table in order to create a custom aggregation (like the Macro Region shown on this picture)

when I go back to the map everything seems to be screwed up, and this is the result I am getting (e.g. islands look like little red crosses, thick polygon borders):

Am I doing something wrong or is it a bug?


texts: Clipping with iterate feature class is only returning some clipped shapefiles with wrong names

texts:
I am trying to create a tool that will clip multiple features (point, line, and polygon) in a geodatabase to a drainage area. I have it set up like this:

Since I have to set this up for someone else to use it, I am testing the tool on a simple geodatabase of national roads, counties, congressional districts, and state capitals (clipping to a state boundary). However, when I run the tool, it only produces two files in the destination geodatabase, _cliproads and _clipstate_caps (the clipped roads and state capitols file). And when I examine the _clipstate_caps layer, it is a polygon layer instead of a point layer, and it contains the clipped features of the congressional districts layer.

Can someone tell me what I am doing wrong?

The model runs fine from the ModelBuilder screen. It does everything it is supposed to. Problems only arise when I try to run it from the dialogue box.


texts: Reaching cm (inch) level elevation precision using drones

texts: A road construction firm needs altitude data of road path before and after the construction and the error should be less than 3 cm (2 inch). Normally, they put a ground station and workers manually read altitude data every 10 meters. How this precision would be achieved using drones? The RTK receiver on drone only tells the exact location of drone, not the location (altitude) of the ground below the drone. One way that came to my mind is combining RTK+LiDAR. The RTK gives altitude of drone and laser tells its distance to ground. However, the error in RTK plus error in laser would more likely reach above 3 cm. Thus, how the drone companies solved this issue?


texts: Band value selection in Google Earth Engine

texts: I am working on Sentinel-5p NO2 data in Google Earth Engine. The band value of 'NO2_column_number_density' starts from negative numbers i.e.(Low : -1.34484e-05) but I need low value to be start from 0. Is there any filter in Google Earth Engine to start band values from 0.


texts: QGIS style polylines with alternate colours

texts:
I have a bunch of shapefiles with multiple polylines. Each consistst of 5 up to 30 polylines. I would like to define a style which colours them alternating two colours.

When I use the Categorized style and define a colour ramp of two colours I get the first 10 polylines in blue and the next 10 in red. Is there an easy way to do this or should I manually create a rule-based style in which all odd numbers are red and all even numbers are blue?


texts:
I have the following gdal_retile command I am trying to run, but it errors on the --optfile:

gdal_retile.py -v -r bilinear -levels 8 -ps 512 512 -co "tiled=YES" -targetDir pyramid --optfile files.txt

   Unrecognised command option: -
   Usage: gdal_retile.py
        [-v] [-co NAME=VALUE]* [-of out_format]
        [-ps pixelWidth pixelHeight]
        [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
               CInt16/CInt32/CFloat32/CFloat64}]
        [ -tileIndex tileIndexName [-tileIndexField fieldName]]
        [ -csv fileName [-csvDelim delimiter]]
        [-s_srs srs_def]  [-pyramidOnly] -levels numberoflevels
        [-r {near/bilinear/cubic/cubicspline/lanczos}]
        [-useDirForEachRow]
        -targetDir TileDirectory input_files

Any idea why this does not work?