texts: OpenLayers: Showing a WFS overlay on top of Google Maps

texts:
I have found similar questions asked previously, but I still fail to get a simple WFS point layer to work using OpenLayers and Google Maps. The original code I used was based on WFS polygon layers I have done previously, I then re-created a page based on the OpenGeo Introduction and other sources, all without success.

The WFS side of things seems fine -- I can see the WFS traffic and it contains co-ordinates in the target space. But I can't see anything rendered at all. Even adding a feature via code doesn't get rendered, in the example code below the feature in the program does display if added to the "vectorLayer", but not if added to the "wfsLayer".

On the risk of appearing stupid: what is wrong with this code?

    <title>OpenLayers: Google Layer Example</title>
    <link rel="stylesheet" href="http://openlayers.org/api/theme/default/style.css" type="text/css" />
    <link rel="stylesheet" href="http://openlayers.org/dev/examples/style.css" type="text/css" />

Note that I replaced some values with "XXX" as to not disclose the client project. As I said: the WFS traffic looks fine, so I assume these values are not relevant for the question. Points in the WFS come through like this:

<gml:Point srsDimension="2" srsName="http://www.opengis.net/gml/srs/epsg.xml#900913">
    <gml:pos>1.2884006544922689E7 -3801543.4695909135</gml:pos>
</gml:Point>

which seems fine for the spherical Mercator projection around the target location.


texts:
When creating a multi-ring buffer I see two methodologies:

Create buffer around feature, buffer the resulting buffer, and repeat for as many rings as required.
Using the initial feature to create the buffer, and use the same initial feature to create subsequent rings, just increasing the buffer distance.

Will these produce the same results? If not which would be the preferred methodology?

This will be implemented in pyQGIS.


texts:
I have downloaded some shapefiles from two different sources both of which show different aspects of the same area.

However, they have different Coordinate reference systems. One has the reference system GCS_Airy 1830, while the other has OSBG_1936_British_National_Grid.

Will this prove to be an issue down the line, and if so how do I choose between the two systems?


texts: Debugging InvalidCastException was unhandled by user code error from ArcObjects DataStat?

texts:
I have some code that adds a bunch of group layers into ArcMap based on a table within ArcMap.  The code works fine in VBA, however in VS 2010 I get an error on the "pEnumVar = pDataStat.UniqueValues" line (see code below).

  ERROR:

  InvalidCastException was unhandled by user code

  Unable to cast object of type
  'System.Runtime.InteropServices.CustomMarshalers.EnumeratorViewOfEnumVariant'
  to type 'ESRI.ArcGIS.esriSystem.IEnumVariantSimple'.

I'm guessing .NET does not like to use the IEmumVariantSimple object?  Does anyone have a suggestion of how I can work around this or update this line with a different enumeration object?

        Dim pEnumVar As IEnumVariantSimple, value As Object

        Dim pODGSLyr As New ODGSLayer

        'Start
        'Find the Layer Files metadata table
        pStdTableColl = m_map
        For i = 0 To pStdTableColl.StandaloneTableCount - 1
            pStdTable = pStdTableColl.StandaloneTable(i)
            If pStdTable.Name = "ODGSLAYERS" Then
                pTable = pStdTable
                lngFldLayerName = pTable.FindField("LAYERNAME")
                lngFldPath = pTable.FindField("PATH")
                lngFldGroupOrder = pTable.FindField("GROUPORDER")
                lngFldGroupName = pTable.FindField("GROUPNAME")
                lngFldGroupTOCOrder = pTable.FindField("GROUPTOCORDER")
                lngFldGroupVis = pTable.FindField("GROUPVISABLE")
            End If
        Next i

        If pStdTable Is Nothing Then
            Exit Sub
        End If

        'Sort the Table
        pTableSort = New TableSort
        With pTableSort
            .Fields = "GROUPTOCORDER, GROUPNAME"
            .Ascending("GROUPTOCORDER") = True
            .Ascending("GROUPNAME") = True
            .QueryFilter = Nothing
            .Table = pTable
        End With
        pTableSort.Sort(Nothing)

        pCursor = pTableSort.Rows

        'Find Unique Values in the Table
        pDataStat = New DataStatistics
        pDataStat.Field = "GROUPNAME"
        pDataStat.Cursor = pCursor

        pEnumVar = pDataStat.UniqueValues '<<<<<<<<<< ERROR HERE
        value = pEnumVar.Next

Update 1

Update 1 is an update from @Kirk Kuykendall comment.  The new error is on the "lngGroupTOCOrder = pRowLyrs.Value(lngFldGroupTOCOrder)" line (see code below).  I'm not sure what I need to do next.

ERROR:

NullReferenceException ws unhandled by user code
Object reference not set to an instance of an object.

 pEnumVar = pDataStat.UniqueValues
        value = pEnumVar.MoveNext
        Do Until IsNothing(value)
            'Now resort the table based upon the layer order in the group
            pQf = New QueryFilter
            pQf.WhereClause = "[GROUPNAME] = '" & value & "'"
            pLyrCursor = pTable.Search(pQf, False)
            pTblSortLyrs = New TableSort
            With pTblSortLyrs
                .Fields = "GROUPORDER"
                .Ascending("GROUPORDER") = True
                .QueryFilter = pQf
                .Table = pTable
            End With
            pTblSortLyrs.Sort(Nothing)

            'Get the newly sorted rows and create the new Group and Layers inside the Group
            pLyrCursor = pTblSortLyrs.Rows
            pRowLyrs = pLyrCursor.NextRow

            'Create the new Group
            lngGroupTOCOrder = pRowLyrs.Value(lngFldGroupTOCOrder) '<<<<<<< NEW ERROR
            blnGroupVis = pRowLyrs.Value(lngFldGroupVis)


texts: Meaning of red cross icon next to Geoprocessing tool/model?

texts:
I am using ArcGIS 10.1 for Desktop.

When I double-click on the toolbox command line, it is not working and appears as a red cross symbol in command tool icon.

So please, what can I do?


texts: Pythonic way to recalculate feature extent in ArcGIS 10.2

texts: I am creating a tool that will automate a tedious process of data standardization. One of the steps is to recalculate the feature extent of a SDE file. Surprisingly I cannot find anything useful in regards to recalculating feature extent, except for a VBA script here: Recalculating full extent of File Geodatabase feature class? and a tool here: http://www.arcgis.com/home/item.html?id=576ebff497094274bdd7c76e456de6c3. Unfortunately I don't know VBA well enough to understand fully what is going on or how to incorporate it into my existing Python script, nor do I know how to call an add-in. Also I will not be the only one using this tool and I'd rather avoid making everyone who uses the tool download the add-in. I've heard about reboxing but it doesn't seem like it's used anymore. Does anyone have any ideas on what I could do or how to use the VBA script in my Python script?


texts: Is it possible to read shapefile field aliases values using fiona?

texts: I am using fiona to read a shapefile, but not sure about how to get field aliases. Is this possible using this great library?


texts: Multi-Otsu algorithm in Google Earth Engine

texts:
I was able to implement a single Otsu threshold using the following code:

// Return the DN that maximizes interclass variance in B5 (in the region).
var otsu = function(histogram) {
var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
var size = means.length().get([0]);
var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
var mean = sum.divide(total);

var indices = ee.List.sequence(1, size);

// Compute between sum of squares, where each mean partitions the data.
var bss = indices.map(function(i) {
var aCounts = counts.slice(0, 0, i);
var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
var aMeans = means.slice(0, 0, i);
var aMean = aMeans.multiply(aCounts)
    .reduce(ee.Reducer.sum(), [0]).get([0])
    .divide(aCount);
var bCount = total.subtract(aCount);
var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
       bCount.multiply(bMean.subtract(mean).pow(2)));
});
print(bss)
print(ui.Chart.array.values(ee.Array(bss), 0, means));

// Return the mean value corresponding to the maximum BSS.
return means.sort(bss).get([-1]);
};

var threshold = otsu(histogram.get('B5_histogram'));
print('threshold', threshold);

var classA = image.select('B5').lt(threshold);

Map.addLayer(classA.mask(classA), {palette: 'blue'}, 'class A');

Now I want to convert this algorithm into Multi-Otsu thresholding algorithm. Any idea on how to do that in Google Earth Engine?


texts: Time Slider feature in ArcGIS 10.5 is displaying layer time extent to 1/1/0001 12:00:00 AM

texts:
I've add a polygon shapefile that has a field that contains date and time in it's column.

However, when I attempt to calculate the layer time extent I'm getting the values of 1/1/0001 12:00:00 AM

Here is a screenshot of what I'm seeing. I've searched around and haven't really found anything specific on this item, like

Time Slider in ArcMap 10 and Time Format Issues;
Setting the time properties on data

I'm sure it is something simple, but I thought I might put this out here to see if anyone may have encountered previously.


texts: Creating a map highlighting areas of potential visibility of a high structure

texts:
The end goal is to add a point which will represent a mast/pylon of x height and for QGIS to generate polygons for the surrounding area which will have direct visibility of the top of the structure from ground level. (I'm using version 2.18.5). Woodland and buildings could make things more complicated, so to simplify we can just stick to the terrain/topography which leads to impeding visibility and tackle woodland at a later stage.

The image below is a perfect example of what I'm attempting to accomplish. It shows a nuclear power plant planned for the west coast of England in Cumbria. The highlighted yellow areas are zones of potential visibility. Is QGIS capable of such a thing? Browsing through the figures released for the power plant, it seemed they were all created in GIS software. The OS mini greyscale map being used as their base map in this case.


texts: Using Stamen base layer offline in mobile application?

texts:
I am creating a mobile application using html/css/js and I am in the process of making it operational offline only.

I am using Stamen as a base layer so I want to use it offline too.

If this is possible, then how?


texts: Coded Value Domain is showing, but not Description?

texts:
I have a Domain that is working, but the issue is when I query on the feature.  I need to be able to search on the Coded Description, not the Code Value.

The attribute table is giving me the Coded Description dropdowns, but when I try to query on any of these values, I need to use the code.

No subtypes on this feature either.


texts: Where can I find a comprehensive overview of QGIS geoalgorithms?

texts:
i am looking for an overview of processing toolbox / qgis geoalgorithms tools
unfortunately the link http://www.qgis.org/de/docs/user_manual/processing/toolbox.html seems to be dead.
SEXTANTE ... (which apparently (?) is the ancestor) doesn't provide any up to date documentation

does anyone have any idea where to find a comprehensive tool-reference including detailed description of in- out-put and parameters ?


texts: Prevent clusters refresh when modifying the style of a feature in OpenLayers?

texts:
I am increasing the size of a feature on mouse hover, but there is a significant delay because, as I've noticed in a browser Performance analysis, all the clusters are refreshed and re-rendered when the feature's "changed" event is trigger.
I've tried scaling the feature in two ways:
marker.getStyle()[0].getImage().setScale(1 * growthRate);
marker.getStyle()[1].getImage().setScale(1.5 * growthRate);
marker.changed()

and
var bg = marker.getStyle()[0];
var icon = marker.getStyle()[1];
bg.getImage().setScale(1 * growthRate);
icon.getImage().setScale(1.5 * growthRate);
marker.setStyle([bg, icon]);

This is how I initialize my layers and the map:
      this.markerVectors = new ol.source.Vector({features: []});

      var clusterSource = new ol.source.Cluster({
        distance: 50,
        source: this.markerVectors
      });

      var clusterCache = {};
      var clusters = new ol.layer.Vector({
        source: clusterSource,
        style: function (feature) {
          var size = feature.get('features').length;
          if (size === 1){
            return feature.get('features')[0].getStyle();
          }
          var color = '3399CC';

          feature.setProperties({ftype: 'cluster');
          var textColor = '#' + utils.invertHex(color);
          color = '#' + color;
          if (clusterCache[size] !== undefined){
            return clusterCache[size];
          }
          var radius = 20;
          if (size < 10) radius = 12;
          if (size >= 10) radius = 15;
          if (size >= 100) radius = 20;
          if (size >= 1000) radius = 30;
          style = new ol.style.Style({
            image: new ol.style.Circle({
              radius: radius,
              stroke: new ol.style.Stroke({
                color: '#000000',
                width: 2
              }),
              fill: new ol.style.Fill({
                color: color,
              }),
            }),
            text: new ol.style.Text({
              text: size.toString(),
              fill: new ol.style.Fill({
                color: textColor//'#fff',
              }),
              scale: [1.2,1.2]
            }),
            zIndex: 8
          });
          clusterCache[size] = style;
          return style;
        },
      });

      this.mapObj = new ol.Map({
        view: new ol.View({
          center: this.center,
          zoom: 8
        }),
        layers: [
          new ol.layer.Tile({source: new ol.source.OSM()}),
          clusters
        ],
        target: 'map'
      });

On smaller input (<100 features) the delay is not noticeable, but on usual input of >1000 features the delay is about 200-500ms. Is there a way to update a feature without triggering the clusters' update function? Or at least to trigger only the cluster in which the feature is present (although it won't modify the state/style of the cluster).


texts:
I successfully merged my imagery data, but the problem is that "merge" function treated an alpha as a transparency instead of an additional channel which has to be merged.

How to merge all rgba properly?

Unfortunately there is no -extend rgba function.


texts: Installing OpenGeo suite breaks tomcat on an instance of AWS Ubuntu Server 12.04

texts:
I am running an instance of Ubuntu Server 12.04 on AWS and after installing the OpenGeo suite have been unable to get access to tomcat6 via localhost:8080 and thus no access to the OpenGeo dashboard.

What I have tried so far: allowed access to port 8080 in AWS, disabled the firewall locally, reinstalled OpenGeo numerous times, and used netstat to see that tomcat is not listening. I have installed Apache and it does listen to port 80.

What is interesting is that I can access tomcat when installed without the OpenGeo suite but once I install the suite I can no longer access it. I am a Linux newbie and at a loss. If necessary, please let me know which logs would be helpful to post. Thanks.

Update: Well no luck after almost two weeks of trying. I have come to the conclusion that OpenGeo/OpenPlans intentionally cripples the Community edition of OpenGeo Suite on Amazon-EC2 in order to sell their packaged for a fee version. Fair enough but I can't afford to pay them and Amazon while I learn.

Update 2: Solved (thank you Don Meltz). Because of the limited memory available in a micro-instance a swap file must be created for tomcat to work. My apologies to OpenGeo. I hope this helps someone else with this problem.


texts: ERDAS Imagine - missing colour column on greyscale image

texts:
I'm having trouble with a greyscale NDVI image that I've imported to EI. I want to set the colours manually but there is no colour column to work with.
When I import the file it just loads directly, no loading window opens and there is no option to load as pseudocolour.
The image was generated from a satellite image via a model that was supplied to me. I've used stretch values of both 2 and 8 (as recommended by my tutor) and it makes no difference.
A note I have from doing this previously says "go to the raster group/table tab/column properties and check the box marked editable next to histogram" but that just makes EI lock up.
Any suggestions?


texts:
I am working with OSM road data for a developing nation. Included in the dataset are various small networks of roads that are orphans or islands i.e. not connected to the main road network.  For example, four islands / orphans are evident in the section of roads below.

My aim here is to have all the roads connected that I can run a Steiner analysis without getting errors from nodes that do not connect.

How can I get QGIS to select these so that I can either delete them or connect them to other sections of road network?


texts: Choosing statistic/process to use in ArcGIS Desktop

texts:
I am trying to determine the best process to analyse some predicted water levels against observed values at specific points.  I can do this simply as 1 to 1 relationship to get a r square.  I would however like to analyse this relationship with a third variable such as the distance from the closest surface water (ie lake).

Is there an ArcGIS Desktop 10.1 statistics tool or process that anyone can recommend?


texts: Discoloration of raster edges on export with QGIS?

texts:
I have a raster tint that looks like this in QGIS:

But after exporting as an image through composer the edges look quite different:

The raster is originally a TIF, and I've tried exporting to JPG and PNG through QGIS composer. I'm guessing that the values are being changed on the export/conversion, but I don't know why.

Can someone please explain what is going on and how I can prevent it?


texts: How to attach lines to points QGIS using Postgis?

texts:
My problem is that I have lines that I want to lengthen and tied them to a point

And in case that I have two line one in continuous and the other is cutted I want only the cutted one to be tied( attached or snaped)  to the point and the other to stay at its place.

I've tried this code but in the result the two lines has snapped to the point.

SELECT f.gid as gid,
   ST_Snap(f.Geom, g.Geom, 2) as geom
FROM cablen as f,
    (SELECT ST_Collect(Geom) as Geom FROM boiten) as g

In the last case the line has snapped to the point and become tow so like you see i have three line and in the result i found four ones but I want them to be like this

Update

after editing my code to this ,i Found the result geom Null

SELECT f.gid AS gid,
   ST_Snap(
       ST_StartPoint(f.geom),
       ST_ClosestPoint(g.Geom,  ST_StartPoint(f.geom) ),
       ST_Distance(
           ST_StartPoint(f.geom)  ,
           ST_ClosestPoint(g.Geom, ST_StartPoint(f.geom) )
        ) * 1.2
    ) AS geom
FROM (
   SELECT ST_Collect(Geom) as Geom
   FROM bt_sout
   ) AS g,
   cablen AS f


texts: How to remove black background from geoserver imagepyramid?

texts:
I have a large ImagePyramid published in Geoserver (v 2.4). Each granule in the pyramid has nodata set to 0. I'm confused as to which layer parameter in GeoServer should I change to get a transparent background. I tried changing BackgroundValues; InputTransparentColor and OutputTransparentColor but nothing seemed to work.

How can I make the background transparent?


texts: Raster intersect and gIntersection give different results in R

texts:
I'm trying to clip a polygon shapefile (species distribution) with another polygon grid. The goal is to calculate the area occupied by the species within each grid cell. I was able to do it with gIntersection from 'rgeos' package, however this function does not retain the arributes from the original shapefiles. I've also tried it using intersect from 'raster' package, and although it does keep all the attributes there is an obvious problem wiht the resulting shapefile.

Here is my code.

# libraries used
library(rgdal)
library(raster)
library(rgeos)

# Read shapefiles
mdau <- readOGR(shp.dir, 'myo_dau') # species disribution data
grid <- readOGR(grid.dir, 'grid300km') # grid

plot(mdau, axes=T); plot(grid, add=T) # check shapefiles are correctly ploted

# Clip using intersect from 'raster' package
inter.raster <- intersect(mdau, grid)
# here I get warnings, which I suspect may be the cause of the unexpected result
# Example Warning messages:
# 1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Ring Self-intersection at or near point -354560.80502695002 8127597.0226868102

inter.raster$area <- area(inter.raster) # area calculation

# Clip using gIntersection from 'rgeos' package
inter.rgeos <- gIntersection(mdau, grid, byid = TRUE)
inter.rgeos$area <- area(inter.rgeos) # area calculation

# Plot results
par(mfrow=c(1,2))
plot(inter.raster, axes=T, main='Raster::Intersect'); plot(inter.rgeos, axes=T, main='Rgeos::gIntersection')

This is the compared output

I would be happy with gIntersection output, however it does not keep attributes from the original shapefiles.

head(inter.rgeos@data)
        area
1   82349136
2   82349136
3   41174568
4   82349136
5 1297047409
6  150000000

Any idea or suggestion of what might be going on?


texts:
I'm having trouble getting a PostGIS raster to load and display in MapServer.  I'm not sure, but suspect it's a version problem.  My mapfile layer and the error message are shown below.

LAYER
    NAME bio12_2_5m
    TYPE raster
    STATUS ON
    DATA "PG:host=localhost port=31600 dbname='wdata' user='postgres' password='*****' schema='worldclim' table='bio12_2_5m'"
    PROCESSING "NODATA=-9999"
    PROCESSING "SCALE=AUTO"
END

msDrawMap(): Image handling error. Failed to draw layer named 'bio12_2_5m'.
msDrawRasterLayerLow(): Unable to access file. Corrupt, empty or missing file 'PG:host=localhost port=31600 dbname='wdata' user='postgres' password='**' schema='worldclim' table='bio12_2_5m'' for layer 'bio12_2_5m'. Sorry, but table bio12_2_5m doesn't seem to have regular blocking                 arrangement. Only tables with regular blocking arrangement can                 be read from now

This is WorldClim data loaded into PostGIS with raster2pgsql like:

raster2pgsql -s 4326 bio12.bil worldclim.bio12_2_5m > bio12.sql

To get around what seems like a blocking problem I've also tried:

raster2pgsql -s 4326 -C -r bio12.bil worldclim.test > test.sql
raster2pgsql -s 4326 -t 10x10 bio12.bil worldclim.test2 > test2.sql

And with the MapServer layer I've tried setting mode=1 or mode=2, still no success.

Versions ... Ubuntu 12.04 LTS, PostgreSQL 9.1.11, GDAL 1.10.0, PostGIS 2.0.4 with raster support, and MapServer 6.0.1.  But (always a but) to be fair I started with GDAL 1.7.3 via the default ubuntu repository, and at the end updated GDAL to 1.10.0 via the Ubuntu GIS "unstable" repository since it seems like the newer versions of GDAL are required to address/support the blocking items.

Thanks, very close I suspect and not sure if installing everything from scratch would help starting with the 1.10.0 GDAL.  Also tried building GDAL from source and no success there.


texts: cropToContents not working properly in PyQGIS

texts:
I have a template qpt which I read via pyqgis script to export multiple maps as png. I have a map and a legend in qpt file and I chose "Only show items inside linked map" so the legend shows the items only in the map. While exporting images, I would like to crop to content so that I do not have extra empty white space. Therefore, I am using image_export_settings.cropToContents = True. However, it seems not working. This is how my map looks when I export

And if I used image_export_settings.cropToContents = True, it looks like. It just crops the legend as well.

Part of my code looks like this
# Read the layout temp
templateFile = open(xxx_qpt, 'r')
templateContent = templateFile.read()
templateFile.close()

doc = QDomDocument()
doc.setContent(templateContent)

layout = QgsLayout(project)
layout.initializeDefaults()

pc = layout.pageCollection()
pc.page(0).setPageSize('A4', QgsLayoutItemPage.Orientation.Portrait )

items, ok = layout.loadFromTemplate(doc, QgsReadWriteContext(), False)

# Adjust frame
mapFrame = layout.itemById("Map1")
mapFrame.zoomToExtent(bbox[0])

# Get the scale of  Layout
referencemap = layout.referenceMap()

wmts_scale = [
    1126,
    2252,
    4504,
    9008,
    18016,
    ]

layer_extent_scale = referencemap.scale()
closest_scale = min(wmts_scale, key=lambda x:abs(x-layer_extent_scale))

index = wmts_scale.index(closest_scale)
mapFrame.setScale(wmts_scale[index+3])

# Export an image of the layout object
exporter = QgsLayoutExporter(layout)
image_export_settings = QgsLayoutExporter.ImageExportSettings()
image_export_settings.cropToContents = True
exporter.exportToImage(imgDir + f"\{str(id)}.jpg", image_export_settings)

What am I missing?


texts: Exporting vector XYZ tiles using pyqgis on QGIS 3.18+

texts:
I am trying to automate the exporting of vector tiles using PyQgis. Without automation, using the tool in the processing toolbox called Write vector tiles(XYZ) under the vector tile category, it creates the output that i need but when clicking on the help button it leads to a 404 page, so that doesnt help.
This is the only code that i have found but it doesnt produce any xyz file directories.
QgsVectorFileWriter.writeAsVectorFormat(osmlayer, outputPath, 'utf-8', osmlayer.crs(),"CSV",layerOptions=['GEOMETRY=AS_XYZ'])

This is my full code
from qgis.core import *
from PyQt5.QtGui import *
import processing
import os
url ='C:/Users/wanza/Desktop/GISOMv1/QGISDATA/qgisdata'
output='C:/Users/wanza/Desktop/GISOMv1/QGISDATA/output/swkosm.shp'
foldlist=os.listdir(url)
url=url+'/'+foldlist[0]
foldlist=os.listdir(url)
for x in range(len(foldlist)):
    bob = foldlist[x]
    if(bob[-3:]=='shp'):
        rob=bob
url=url+'/'+rob
layer = QgsVectorLayer(url,'temp','ogr')
QgsProject.instance().addMapLayer(layer)
ext=layer.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
coords = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax)
osmdata = 'C:/Users/wanza/Desktop/GISOMv1/QGISDATA/qgisdata/3qosmdata/gis_osm_buildings_a_free_1.shp'
processing.run("gdal:clipvectorbyextent",{'INPUT':osmdata,'EXTENT':coords,'OPTIONS':'','OUTPUT':output})
osmlayer = QgsVectorLayer(output,'temposm','ogr')
QgsProject.instance().addMapLayer(osmlayer)
outDir = 'C:/Users/wanza/Desktop/GISOMv1/QGISDATA/swkosm'
if not os.path.isdir(outDir):
    os.mkdir(outDir)
outputPath = outDir
QgsVectorFileWriter.writeAsVectorFormat(osmlayer, outputPath, 'utf-8', osmlayer.crs(),"CSV",layerOptions=['GEOMETRY=AS_XYZ'])

Mind the code at the top as it is mostly hard coded just to confirm the functionality of what I am doing so far. The code to export from what i have found is at the bottom, and it currently doesnt work. Just assume that i need a vector layer to export into xyz tiles. I just cant find the code for it.


texts:
I have a simple tasks that regarding combining multiple Postgis tables. Let's say I have 4 tables called tx_shape, ca_cd_1_shape,wy_shape, and ma_suffulk_co_shape. Two of which are states, one is a congressional district, and another is a county.

I want to combine them into one giant shape table. I tried using referring to this post, and ran

CREATE TABLE merged AS (
SELECT gid, geom FROM tx_shape
UNION
SELECT gid, geom from wy_shape
UNION
SELECT gid, geom FROM ca_cd_1_shape
UNION
SELECT gid, geom FROM ma_suffolk_co_shape);

But all it returns is a blank. The gid are displaying correctly but when I export it using pgsql2shp with

pgsql2shp -f "test_output" -h localhost -u username -P password my_db "select * from merged"

All I get is a blank, am I missing something here? It should be relatively straightforward.


texts: How to skip some steps, and continue the running of a model from a later point, bases on a condition?

texts:
I am building a Model in ArcGIS Desktop 10.8.
At some point in my model, I am dealing with Line feature classes. I turn them into Polygon feature classes, and then I merge them into one feature class.
The source of my problem is, that some of the inputs might be empty. If I run the Feature To Polygon tool with these inputs, it gives me the following error message:
Error executing function. The geometry is not M-aware Failed to execute
I tried to solve the problem, based on this thread: Geometry is not M-Aware?. It didn't help. I am working in a geodatabase, and the M-aware option is turned off in the environmental settings.
I thought that preconditioning would work, so I tried (as you can see in the 3rd and 4th row of the attached image). The problem here is that it won't run the Feature To Polygon tool indeed, but it won't run anything after this tool either. So it won't eventually merge the other, not empty feature classes either.
Is there a way to skip some steps, based on a condition, and continue the running from an other point?
In my case, can I skip the Feature To Polygon and the Calculate Field tools if my feature class is empty, but then run the Merge tool with the other, not empty feature classes?


texts: What does it mean that shapefiles are "mostly" open specification?

texts:
Esri Shapefiles are commonly considered the standard format for GIS vector data.

I noticed on Wikipedia and the OpenStreetMap wiki it is noted (emphasis mine):

  The shapefile format is a popular geospatial vector data format for
  geographic information system (GIS) software. It is developed and
  regulated by Esri as a (mostly) open specification for data
  interoperability among Esri and other GIS software products.

The citation for this is the ESRI shapefile technical description, which doesn't make mention of the term "open specification".

The Wikipedia definition for an open specification is:

  An open specification is a specification created and controlled, in an
  open and fair process, by an association or a standardization body
  intending to achieve interoperability and interchangeability.

It also cautions against confusing open standards and open specifications, which perhaps I'm doing.  Is the shapefile a standard or a specification?  Does the Esri shapefile meet the above definition, and if not, what does it mean to say the shapefile is an open specification?  Moreover, what does it mean that it is "mostly" open?  What are the limitations to this specification's openness?


texts: How can I find the elevation a set distance upstream and downstream from a layer of crossing locations

texts:
I have a channel network which I've intersected with a proposed road network to produce crossing locations. I need to find the channel slope averaged over a 20m stream length.

I also have a raster of elevations

Is there a way to 'look' upstream and downstream the channel network from a set of crossing locations?


texts:
I have arcgis 10.3 and I am trying to run a task using python but the python window within arcmap won't import arcpy.  This never happened to me with 10.2.2.  The error returned is:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 24, in <module>
    from arcpy.toolbox import *
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\toolbox.py", line 15, in <module>
    from .ddd import ASCII3DToFeatureClass as ASCII3DToFeatureClass_3d
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\ddd.py", line 23, in <module>
    from arcpy.arcobjects.arcobjectconversion import convertArcObjectToPythonObject
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\__init__.py", line 19, in <module>
    arcobjectconversion.initialize_conversion_mapping()
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\arcobjectconversion.py", line 26, in initialize_conversion_mapping
    from .._mapping import MapDocument, PageLayout, DataFrame, DataFrameTime, DataDrivenPages, \
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\_mapping.py", line 19, in <module>
    from arcpy.utils import logcall, ArgAdaptor
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 17, in <module>
    import inspect
ImportError: No module named inspect


texts:
I am running QGIS 2.18.7 on Mac OS (10.11.6).
The application crashes frequently. I usually only manage a few minutes work (while saving all the time) before it crashes.
I recently upgraded to this version as the last version was crashing too often... glad that helped!
There doesn't seem to be any common thing that causes this. Sometimes it crashes when I try to edit a shapefile, sometimes when I open a raster layer, sometimes when I zoom.... endless number of things.
I have tried setting up a new workspace, and keeping it trim and none of my open files are exceptionally large.
I have tried uninstalling all the plugins.

I have also checked other threads and noticed a few people having these issues, and the only thing that helped anybody was to follow the advice on the install page to rename or move the location of 'Qt4 developer components'. However I have no idea what these are.. or where to find them!! I search in finder and look in the usual places but no sign of these items.


texts: Creating a weekly FGDB backup of ArcSDE?

texts:
I need to create a script that will copy our county's parcel data, maintained in SDE/SQL, to a file geodatabase.  This script will be run as a scheduled task over the weekend.  What GP tool do you recommend for this process? Consider that our schema is 'set in stone' and we have no relationship classes to maintain.  I have tested both CopyFeatures_management and FeatureClassToFeatureClass_conversion.  Both work.  FC to FC is slightly faster, but speed is not really an issue.  In theory, CopyFeatures, run with Overwrite = True, matches our workflow as we are not converting to a shapefile or coverage.  FeatureClassToGeodatabase_conversion may be the best choice, as there are about 25 feature classes that need to be copied to the same FGDB.

Also, I need to ensure that the machine running the script is able to connect to the SDE data.  We are using DBO-schema with OS authentication. The script will likely be run by the same server that houses our SDE.  Does anyone have experience with a similar scenario?


texts: How to find the "first" intersection of a line with anything else

texts:
I am trying to find the first intersection of a line with anything else.

By first intersection I mean the first intersection, if you follow the line in its direction from start to end.

Check the image below, it should make it clear. Black is a line, blue is another geometry, here a polygon, red is the point I am looking for.

Maybe there is some standard algorithm for this that I am not aware of? Or maybe a combination of multiple ones?

I am not interested in solutions of closed-source software.


texts: In QGIS symbology controls there is a "Point Displacment" option for offsetting points that are close visual clarity.  Is there a similar control for MapServer?


texts: Form opening after stopEditing pyqgs Qgis

texts:
I made a form that binds some fields on a shp in Qgis (works fine), and after stop editing event the form appears again..
Do you have any ideas?
Made so far:

## action
self.editingObject = QAction(QIcon(":/plugins/Count/icons/pen.png"), "&Editing object", self.iface.mainWindow(), toolTip ='<b>Editing object</b>\n Editing object',)
self.editingObject.triggered.connect(self.editing)

def editing(self):

    if not iface.activeLayer():
        QMessageBox.about(None, 'Atention', 'No layer selected')

    self.layer_activ = iface.activeLayer()
    LayerSource = self.layer_activ.source()
    sufix = 'judete_ro.shp'
    if LayerSource.endswith(sufix) :

         self.iface.activeLayer().startEditing ()

         QMessageBox.about(None, 'Info', 'Editing started on layer' + str(iface.activeLayer().name()))

         self.iface.actionAddFeature().trigger()

     if self.editingObject.triggered:
         self.layer_activ.featureAdded.connect(self.select_id)
def select_id(self, Fid):
       self.iface.activeLayer().setSelectedFeatures([Fid])
       self.editfrm = editForm()
       self.editfrm.show()
class editForm(QtGui.QDialog):
    def __init__(self):
       QtGui.QDialog.__init__(self)
       self.ui = Ui_EditForm()
       self.ui.setupUi(self)
       self.ui.pushButton.clicked.connect(self.accept)

    def accept(self):

      for i in iface.activeLayer().selectedFeatures():

         fid = int(i.id())
         provider = iface.activeLayer().dataProvider()
         name = iface.activeLayer().fieldNameIndex("NAME")
         siruta = iface.activeLayer().fieldNameIndex("SIRUTA")
         textNAME = self.ui.NAME.text()
         textSIRUTA = self.ui.SIRUTA.text()
         attr = {name: textNAME, siruta: textSIRUTA}
         iface.activeLayer().changeAttributeValue(fid,name,textNAME)
         iface.activeLayer().changeAttributeValue(fid,siruta,textSIRUTA)
         iface.activeLayer().setSelectedFeatures([])
      self.close()


texts: Updating production ArcSDE database from staging?

texts:
We have development, staging and production instances of the same ArcSDE database, similar to this ancient thread from the ArcGIS Discussion Forums.

Our development to staging deployment employs lots of schema changes with exclusive locks that involve dropping tables and recreating them. Our schema itself is pretty fluid as we are always adding new data to it.

Our staging database is an exact copy of what we want our production database to be.

How do I deploy our staging ArcSDE instance to production while minimizing or eliminating downtime?

We use SQL Server 2012.


texts:
In QGIS, I have a WMS layer. I am saving it to my hard drive by right click > "Save As...". Just like in HDunn's answer here. This is about my understanding of the topic, already this code is strange and not really understandable for me.

However, "Save As..." is resampling the original data by averageing it, which is not OK for me since the original raster is a categorical one. I do not know the resolution of the original WMS file (it is not included in it's description).

How could I save it with either the original resolution (and CRS) and/or with other resampling method than averaging? Natural Neighbour would be fair enough although Most Frequent Value would be most accurate. Is that possible without Python? I also understand of R, if that helps.

The WMS URL is: https://maps.rissac.hu:6443/arcgis/services/dosoremi_web_mercator/MapServer/WMSServer?LAYERS=4&PROJECTION=EPSG%3A3857&TRANSPARENT=TRUE&VERSION=1.1.1


texts: Zooming into center of map on load of ArcGIS API for JavaScript 4?

texts:
I was looking through the ArcGIS code samples and I found this. I really like the way that it zooms in from afar, but I can't seem to figure out how to fly in to the location on load like this version does.

Can anyone help me figure out what lines drive this feature?


texts: Interpolation grid for scattered data having (latitude, longitude) coordinates

texts:
I've got some scattered data in the form of (latitude, longitude, someParameterValue). I'm using inverse distance weighting interpolation method to interpolate them in a rectangular grid of pixels. Presently I'm generating the query points for that grid, in python, as given below. Please note that I've converted the (latitude, longitude) coordinates to cartesian (x, y) coordinates :

        xr = int(math.ceil(xmax-xmin));
        yr = int(math.ceil(ymax-ymin))
        xres=yres=1
        xr = math.ceil(xr/xres)+1
        yr = math.ceil(yr/yres)+1
        npts = int(xr * yr)
        queryPts = np.zeros(shape=(npts,2))
        x1=xmin; y1=ymin
        idx=0
        while(x1 <= xmax):
            while(y1 <= ymax):
                #the 2D array of query points is populated here
                queryPts[idx] = [x1, y1]
                idx += 1
                y1 += yres
            y1=ymin
            x1 += xres

where xmin, ymin, xmax, ymax are the minimum and maximum values of x and y coordinates respectively. Here I feel that populating the query points at intervals of 1 in each of x and y axes is not the right way to go. After calculating the grid of interpolated values, I'm using gdal to turn it into a raster image with the interpolated values scaled to 0-255 for the pixels. A sample image, that I assume to be inappropriate, is shown here.

Would like to get suggestions on the following:

What would be the proper way to generate query points for an interpolation grid i.e. how to set the resolution of the points on x and y axes whose interpolated values that we are going to calculate?
I can see that the above question is related to the pixel resolution that we want in the final interpolated image. Hence the above question could be asked as: how to generate 2D arrays of (x,y) points (query points?) that correspond to pixels in the rendered image?


texts: How to capture contour lines from OSM cycle map as vector?

texts:
In QGIS, there is a plugin to insert a OSM cycle map layer, that brings as an image from http. But I can't convert this OSM Plugin layer to a raster image.
I can't even imagine how to get the contour lines that are showing by the cycle map as a vector.

Is this possible?

If is not possible through OSM, there is another open resource to get this?


texts:
i am trying to compute the time series graph of the particular location i m getting this error no idea what to do i m using "MODIS/006/MOD11A"

var lst = lstImg.select('LST_Day_1km')
             .filter(ee.Filter.date('2000-01-01', '2019-12-31'))
             .map(function(image){return image.clip(table)})
            .sort('CLOUD COVERAGE',false)
             .filterBounds(table);
    var lstVis = {
      min: 13000.0,
      max: 16500.0,
      palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
      ],
    };
    //Map.setCenter(6.746, 46.529, 2);
    Map.addLayer(lst,lstVis,'LST');
    print(ui.Chart.image.series(lst.select('LST_Day_1km'), table, ee.Reducer.mean(),30));

ERROR

Error generating chart: Collection query aborted after accumulating over 5000 elements.

Collection query aborted after accumulating over 5000 elements.

Script link

https://code.earthengine.google.com/ba74a316399d35095a3d8c5a9ccabb10


texts: Openlayers 3, jquery Ajax initialising in Jquery.success()

texts:
I have this nice fiddle: http://jsfiddle.net/mindsupport/cujfjby4/14/

Openlayers 3 map initialized in ajax request and all data (points, line etc.) attached in jquery success() method, but map doesn't seen, it's is not initialized, but all objects a added to the global linkMap variable. Why? I've tried withoiut Ajax call, gust pure json data and it works...


texts: How do I add an offset to a lat/lon to get a new point using PostGIS

texts: Using PostGIS I have a lat/lon point and want to add an east/west and north/south offset in metres to get a new point. Is that possible?


texts: Meaning of ERROR 000539: Error running expression: rcexec() : No module named os?

texts:
In Raster Calculator using:

(FocalSt_meanDEMtif5-FocalSt_minDEMtif1) / (FocalSt_maxDEMtif2-FocalSt_minDEMtif1)

What is the meaning of

  ERROR 000539: Error running expression: rcexec() : No module named os


texts: I'm trying to compute zonal statistics but get an error saying that "negative dimensions are not allowed". I'm assuming its because my raster file extends into negative coordinates. Extent of the file is xmin: 0, xmax: 21600, ymin: -8700, ymax: 0. (The top left corner of the image is 0, 0). I want to move this up so that ymin: 0 and ymax: 8700. Alternatively, if there is another way to avoid this error that is welcome.


texts: Confused again.  I was finally able to run ArcCatalog and create a join in my map document (I will do the deed programmatically eventually) but I was interested in how a join looked once it was created. However, when I iterate thru the layers and read the fields associated to each layer using get_Field(index) in C# I see no joined fields. How can I read (and use) an existing join?


texts: Create Random Points with based on Attributes values in FME

texts:
I have 6 Polygons in one shapefile with attributes Block A to Block F.

Also, I have one Excel file for materials location like shown below screenshot.

My aim is to create a random point within the polygon as per the values mentioned in the table.

For Example:

121 Sand Random points for Block A polygon

6182 Sand Random points for Block B polygon

1 Cement Random points for Block A polygon

210 Cement Random points for Block B polygon

All points if possible one single shp file.


texts: Extract/copy a subset of rasters from a mosaic dataset

texts:
I have a mosaic dataset of 421 images and I want to extract just a few of them to give to a client. Rather than click through the folder selecting the files (because they are not contiguous) I want to transfer is there a way to use the list of selected footprints (using selection features in ArcGIS) to automate the copying? If there is I'm running into a total mind blank and have not found anything close to what I'm trying to do.

Possibly something I could use in Total Commander 8 or CMD prompt (i.e. xcopy/robocopy, etc) because ArcGIS would be too slow I think.

EDIT.

I also have Global Mapper 14 to use of the explorer options aren't feasible.


texts:
I am running this code

import geopandas
df = geopandas.read_file('test.geojson')
ax = df.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')

and I get nothing besides this error

/usr/local/lib/python3.6/dist-packages/matplotlib/backends/backend_gtk3.py:197: Warning: Source ID 8 was not found when attempting to remove it
  GLib.source_remove(self._idle_draw_id)

The GeoJSON is okay, I tested it in QGIS. The same error turns up if I try it with other geofiles.

What am I missing?


texts: How to rename ogr2ogr tables on export to Spatialite?

texts:
I'm exporting data from PostGIS to Spatialite using ogr2ogr. The export is working properly, except I need to be able to rename the tables that are being exported. Here's the command I'm using:

ogr2ogr -f "SQLite" testSQLite.sqlite -progress PG:"host=127.0.0.1 dbname=mydb user=admin password=myPassword active_schema=my_schema schemas=my_schema tables=building_shared,alarm_shared" -lco LAUNDER=yes -dsco SPATIALITE=yes -lco SPATIAL_INDEX=yes -overwrite

This exports a Spatialite database with two tables - building_shared and alarm_shared. I want the two tables to be named building and alarm. I tried adding the following SQL query to the end of the command, with no luck:

-sql "SELECT * FROM building_shared AS building"

How do I rename the tables on export?


texts: Clear resource cache via REST in Geoserver

texts: Is there a way to clear the Geoserver resource cache using REST? Same as "Server status" > "resource cache" > "clear" in the Admin UI.


texts: Geonode to harvest metadata from Geonetwork-like services?

texts: I am in charge of setting up a gis web portal and one of the functionalities would be to harvest some metadata from metadata servers like Geonetwork. I have been thinking of using Geonode because I have already used Geoserver. Apparently Geonode uses Pycsw to publish metadata but I can't see whether it allows to collect metadata from other servers and how? Is it possible to do it through the Geonode interface?


texts: R terra and `Error in x@ptr$geometry(nlyrs) : could not find valid method`

texts:
i am using terra::raster 0.9-1 and, specifically terra::values and the xyz function in terra::rast, but i keep encountering the following
Error in x@ptr$geometry(nlyrs) : could not find valid method
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)
other attached packages:
 [1] terra_0.9-1       gdalUtils_2.0.3.2 geosphere_1.5-10  data.table_1.13.0 ggplot2_3.3.2     dplyr_1.0.2
 [9] rgdal_1.5-16      rgeos_0.5-3       stringr_1.4.0     raster_3.3-13     sp_1.4-4          devtools_2.3.1

# dummy XYZ data
m_xyz <- as.matrix(expand.grid(seq(0.5,999.5,1),seq(0.5,999.5,1)))
m_xyz <- cbind(m_xyz, sample(1:10,1000^2,replace=T))

# coerce to rast (used to be raster::rasterfromXYZ)
rast(m_xyz,type="xyz")
Error in x@ptr$geometry(nlyrs) : could not find valid method

ok, ignore that and just create a rast and try to populate
rBNG <- rast(xmin = 0, xmax = 1000, ymin = 0, ymax = 1000, resolution = 1, crs = "epsg:27700", vals = NA)

# add some values
values(rBNG) <- sample(1:10,ncell(rBNG), replace = T)
Error in x@ptr$geometry(nlyrs) : could not find valid method

rBNG[] <- sample(1:10,ncell(rBNG), replace = T)
Error in x@ptr$geometry(nlyrs) : could not find valid method

i am new to terra, i have moved to sf and stars and now terra is the next step, perhaps there's just something simple i am not understanding
Loading required package: sp
rgdal: version: 1.5-16, (SVN revision 1050)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
Path to GDAL shared files: C:/Users/xxxx/Documents/R/win-library/4.0/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: C:/Users/xxxxx/Documents/R/win-library/4.0/rgdal/proj
Linking to sp version:1.4-2
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

Attaching package: ‘rgdal’

The following object is masked from ‘package:terra’:

    project

Warning message:
package ‘sp’ was built under R version 4.0.3


texts: Why would Nearest Neighbor Analysis give 0 observed mean distance?

texts:
I am working with small datasets (around 100 points or less) with Nearest Neighbor Analysis in QGIS 2.8 (although I also tested in 2.6 and the result is the same). I have the layers referenced with a projected coordinate system.

When giving the results, I always get 0 as Observed Mean Distance, while Expected Mean Distance seems correct. I copy a typical output below:

Observed mean distance:0.0
Expected mean distance:367.212366515
Nearest neighbour index:0.0
N:98
Z-Score:-18.9384277177

Tested the same shapefile in other software and the value is valid (not 0).
Is there any known limitation for the input shapefile? Any workaround?


texts: Calculate coordinates of centroid in arcgis using python

texts:
I have polygon features and want to calculate coordinates of it's centroids in arcgis using python - all goes fine, but it gives me coordinates in meters and not in degrees as i want.

Here is the code:

arcpy.CalculateField_management(folder+"\\"+f+"_fakela_shp.shp", "X", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
arcpy.CalculateField_management(folder+"\\"+f+"_fakela_shp.shp", "Y", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[1]", "PYTHON")

@DECIMALDEGREES should done the job, but it is broken somehow.


texts: Changing CRS in QGIS due to pop up window: "CRS warning! Warning: Input layers have non-matching CRS. This may cause unexpected results.”

texts:
I'm using GIS to map the movement of animals and their speed and density in different habitat types using shape polygons.

When I got to the step ;Count points in polygon' a pop up window appeared saying:"CRS warning! Warning: Input layers have non-matching CRS. This may cause unexpected results.”
I clicked ‘ok' however this meant that the attributes table that I then opened had no GPS speed values.

I'd like to know how to convert the CRS for this habitat shape polygon from 'Selected CRS (EPSG: 4326, WGS 84)' to 'Selected CRS (EPSG: 7405, OSGB 1936 / British National Grid)'.


texts: Clipping raster with an irregular shape?

texts:
How can I clip my DEM raster in a irregular shape in QGIS?

I have used the clipping tool but it only does rectangles I believe.


texts: Dragging feature- linked annotation results in strange snapping and rotation behavior

texts:
I created feature linked annotation for a line feature class. Both the feature and the annotation are feature classes stored within the same feature dataset on an SQL Server Geodatabase.

The initial labeling made use of a simple SQL query which may or may not be relevant. The annotation was created successfully but I'm having a hard time maintaining it due to the following:

When dragging these annotation boxes to fine tune my print product, I use the Edit annotation tool. I've attempted this in both data and layout view (with a focused data frame). When stretching this bubble style annotation it is automatically rotated in a variety of directions and appears to snap to some unknown grid even though all snapping is turned off.  It makes placing annotation impossible.

I posted Feature linked annotation not added when new features are drawn which is a related but different question here today. These are not duplicates.


texts: abind Error when converting a .tif raster data object to an sf object using st_as_sf() in R

texts:
When I try to convert a .tif raster data object to an sf object using st_as_sf() in R, I get the following error:
Error in abind(x, along = length(dim(x)) + 1) : along must be between 0 and 2

Does anyone know how to solve this?
The .tif data for which to reproduce this I got from https://sedac.ciesin.columbia.edu/data/set/gpw-v4-data-quality-indicators-rev11/data-download.
The error occurs both for the Data Context and Mean Administrative Unit Area data. Choose GeoTiff and 30 sec (1 km) resolution.
The reprex R code is then (for the Data Context data, for example):
library(stars)

# Load data
gpw_raw <- read_stars("gpw_v4_data_quality_indicators_rev11_context_30_sec.tif")

# Attemt to convert .tif into sf object
gpw <- st_as_sf(gpw_raw)

One resource I used, which in the end did not allow me to solve the issue, is: https://r-spatial.github.io/stars/articles/stars5.html#vectorizing-a-raster-object-to-an-sf-object
Thanks for your help!


texts: How can i add sub menu under custom top menu in QGIS using python

texts: I have added "My Tools" Menu on top. Now i am clicking on "Plugin" then selecting QGIS Plugin bulider. Here in the menu field a drop down list is coming. How can i add "My tools" in this drop down list.


texts: How to speed up batch watershed delineation?

texts:
I have a large number of pour points - approximately 250,000 for which I want to delineate the upstream watershed area using ArcHydro's "Batch Watershed Delineation". Given that it takes roughly 4-5 hours for 1000 points to delineate on my system I need to find a faster way than just running these as individual processes on my computer. Currently I am testing out running 4 sets of 1000 side-by-side on my 64bit windows i7 system, so that will likely speed up the process, but I think it could be sped up even more through an opportunity to access a "condor" type system which runs multiple applications of a program at any given time.

The catch is that I need to talk to ArcGIS from a command line and also get ArcHydro to work from command line to do the upstream watershed delineation. I know the files that I require, so would likely just need to know if it is possible to extract the code that ArcHydro uses to define the entire upstream catchment area of points. I can't seem to find that easily within the ArcGIS GUI.

Would anyone know of a way to identify/extract the code that ArcHydro is using to run "Batch Watershed Delineation", and if so would you also know how to call to ArcGIS to do that process without opening ArcGIS and running it from a command line?


texts: Where's the GDALTools plugin for QGIS 2.4? It's not listed in list of available plugins, and it's seemingly not installed out of the box. Is it just that it's not yet compatible?


texts: Why does GeoServer WMS not work in OpenLayers? Problem viewing GeoServer layers

texts:
Adding two layers using OpenLayers, one map layer from OSM, second layer adding from GeoServer localhost. But when I load single OSM layer it's showing by default map.

It's okay, but when adding second layer it's not working. it's also disappear old layer map also
and when I remove first layer it's won't show anything from second layer

var MapLayer = new ol.layer.Tile({
            source: new ol.source.OSM()
        });
        var myView = new ol.View({
            center: [74.1296577, 32.2739606],
            zoom: 4,
            projection: projection
        });
        var mymap = new ol.source.TileWMS({
            url: 'http://localhost:8080/geoserver/DHAGIS/wms',
            params: { 'LAYERS': 'DHAGIS:baddoki', 'TILED': true, 'STYLES': 'graves' },
            serverType: 'geoserver',

        });
        var mylayer = new ol.layer.Tile({
            source=mymap
        });
        //////////////////////////
        var map = new ol.Map({
            target: 'map',
            layers: [MapLayer, mylayer
            ],
            view: myView
        });


texts: Query and export a grouped subset of fields from a shapefile layer

texts:
I have several shapefile vector layers that I need to query so I can group records and export a subset of fields to go into an excel spreadsheet.  I thought I could do this by using the graphical modeler and simply using the 'Execute SQL' algorithm on the input layer and then copying the records into excel.  When I try to do this, I get an error that it cannot find the virtual layer ("virtual: Cannot find layer [generated temp file name]").

Can this be done on shapefile layers or do they need to be in a database layer?  Is there a better way to do this?


texts:
Is it possible to draw an arc(curve) in QGIS using Python if  the following parameters are known?

start point of the arc
end point  of the arc
arc length
direction of the arc from start point to end point
Radius (of the circle created from the arc)


texts: Is this possible with ArcGIS using replication?

texts:
I work for a telecom company dealing with fiber management.

We have a sister company south of us who we share some fiber and fiber equipment with.

Both of us will be using ArcSDE databases.

We want to be able to share this area with each other, but we don't want to share our entire system. We want to be able to make changes in this area where we share fiber/equipment and sync the changes back and forth, but only for these areas.

Please see this image to have a visual:

We only want to sync changes made in the red area.

The problem I'm having with using something such as the two-way replica is that if I save the replica into the database where the data currently exists it doesn't merge with the current data, it creates new data with an _1, like instead of appending to the points feature class it makes a new one called points_1. This is a problem because it also does it for the geometric network which is being used by our fiber management software to track the fiber.

Is there anyway to get the replicated data into our original data and still only be able to sync that data?

I'm not sure if this made sense. It's kind of hard to explain.


texts:
I downloaded some shapefiles from the Australian Bureau of Statistics (ABS).
Its basically showing a bunch of suburbs when I view the shapefile through a GIS tool.
The GIS tool that I am using at present is the GIS GRASS 6.4.2 tool.
When I extract the polygons of each postal out of the shapefile I noticed that some postal codes overlap with one another.
I assume this should not be the case, and the correct display is to make sure the overlapping postal codes exclude on another in areas that they overlap.
Is there some tool I can use to make sure that no postal codes overlap one another?


texts: Get "approximate" lat/long coordinate for UTM grid zone

texts: I have data assigned with a UTM grid zone designators, such as 10S for example. Now I need an approximate latitude/longitude pair, preferably the center point of the given grid zone. Can you give me some advice, on how to map those values?


texts: Using update cursor in a while loop to perform analysis on individual polygons gets stuck on first row

texts:
I am working through an open courseware project from PSU found at https://www.e-education.psu.edu/geog485/node/145

I am trying to use an update cursor in a  while loop to perform an analysis on individual polygons within a multipart feature. However, the output is for the first row ad infinitum.

Here is what i have so far:

import arcpy
arcpy.env.workspace = "c:/users/documents/arcgis/arcgis tutorials/psu tutorials/project3/policedata.gdb"
arcpy.env.overwriteOutput = True

# Declare variables for files and fields used in analysis

incidents = "c:/users/documents/arcgis/arcgis tutorials/psu tutorials/project3/policedata.gdb/GraffitiIncidents"
patrol_zones = "c:/users/documents/arcgis/arcgis tutorials/psu tutorials/project3/policedata.gdb/PatrolZones"
incidents_layer = "IncidentsLayer"
patrol_field = "NAME"

cursor = arcpy.UpdateCursor(patrol_zones)
zone = cursor.next()

# Loop through each patrol zone

while zone:
    try:
        arcpy.MakeFeatureLayer_management(incidents,"IncidentsLayer" )
        zone_name = zone.getValue(patrol_field)
        print zone_name

        sql = patrol_field + " = " + "'" + zone_name + "'"

    # Make a feature layer for each zone
        arcpy.MakeFeatureLayer_management(patrol_zones, "CurrentZone", sql)

    # Make selection
        arcpy.SelectLayerByLocation_management(incidents_layer, "CONTAINED_BY", "CurrentZone")
        num_incidents = arcpy.GetCount_management(incidents_layer)

        print "Graffitti incidents = " + str(num_incidents)
##        arcpy.Delete_management("CurrentZone")
##        zone
    except:
        print "Could not perform analysis."

del cursor

What in my code is causing the script to produce results only for the first row?


texts: How do you get the path of a file you saved in memory: or as temp file with pyqgis

texts:
So whenever I do something like buffer in QGIS with Vecotr>>Geoprocessing Tools>>Buffer... the source output I get for my temporary file looks like this:

That works fine if I'm doing more QGIS actions, but then I need to use r.patch and for it to work I need to put actual paths as inputs(or if I don't tell me how to make it accept the source above).

If I do it manually, it knows the path and if I go in the processing history I can see he used:
C:/Users/Phil/AppData/Local/Temp/processing_534fdc910b13414bbab4810d36304212/15b103f620a74a189c2370a8abee2ba0/output.tif

But I have no idea how to go get that path. Right now this is how my buffer looks:

processing.runAndLoadResults("native:buffer", \
{'INPUT':tmp_riv_ponc_Intersect_Clipped,\
'DISTANCE':1.5,\
'SEGMENTS':5,\
'END_CAP_STYLE':0,\
'JOIN_STYLE':0,\
'MITER_LIMIT':2,\
'DISSOLVE':False,\
'OUTPUT':'TEMPORARY_OUTPUT'})

QgsProject.instance().mapLayersByName('Buffered')[0].setName('Riv_Ponc_Intersect_Buffered')
tmp_riv_ponc_Intersect_Buffered = QgsProject.instance().mapLayersByName('Riv_Ponc_Intersect_Buffered')[0].source()

And the tmp_riv_ponc_Intersect_Buffered gives the source of the picture shown above.
I also tried to use the same line but instead of source I put ".dataProvider().dataSourceUri()" but it gave me pretty much the same thing except a couple of lines from the Source were missing.

I also tried to put "memory:" instead of TEMPORARY_OUTPUT but didn't change anything.
I tried to do

output_path = processing.runAndLoadResults(...
{'OUTPUT': 'output_22c6c1b4_1c03_44c8_92dd_3dffc7a6a0d8'} #Output = "memory:"
{'OUTPUT': 'Clipped_5b605e64_620d_4ac4_a1b1_64a329f6c2b5'} #Output = "TEMPORARY_OUTPUT"

They were done on a different processing algorithm, but "memory:" always gives output_... and the temporary gives me the name of the action like Clipped_...

I have found no way of getting the complete path to the temp folder where the files should be.

Thank you!


texts: Why are random characters (e, ?, `, ... ) appearing in my text field (attribute table) of the shapefile in QGIS 3.0?

texts:
I have a shapefile that I have created which has a few features (about 10) and I have input attribute information for each feature and each attribute (set to text).

For example, one of my fields is called Acer negundo (no accents or other odd fonts). Sometimes, the unicode replacement character � appears beside the information in the fields with attribute information and also in the fields with no attribute information.

It is further frustrating because I have a filter (is not NULL) to only show features with attribute information, and the filter doesn't work when the random characters appear.

I've tried changing the coding from system to UTF-8 without any helpful changes.


texts: arcpy.SearchCursor loops one time but then does not repeat

texts:
The part of the script I'm having trouble with is the arcpy.SearchCursor. The cursor runs 1x and then does not re-run. The goal was to re-retrieve the field value and apply it to the image layout. It works on the first run but then applys the same value to every image layout after. The SQL query I'm using limits the records to 1 record, but I was hoping it would re-run and apply the new SQL query.

x = 1

while x <= formatted_rand_ct:
    row_rand = random.randint(0,record_count)
    sqlExp = new_name + "=" + str(row_rand)
    rows = arcpy.SearchCursor(newlayer, sqlExp)
    for row in rows:
        l = row.getValue(field_label)
    arcpy.SelectLayerByAttribute_management(newlayer, "NEW_SELECTION", sqlExp)
    df.zoomToSelectedFeatures()
    if desc.shapeType == "Point":
        df.scale = scale
    if desc.shapeType == "Point" or desc.shapeType == "Polygon":
        elm = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT")[0]
        elm.text = str(field) + " " + str(row_rand)
        elm2 = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT")[3]
        elm2.text = str(l)
        arcpy.mapping.ExportToJPEG(mxd, output + str(row_rand))
        arcpy.RefreshActiveView()
        arcpy.AddMessage("Finished Export of Image: OBJECTID" + " " + str(row_rand))
        x += 1
    else:
        arcpy.AddMessage("Points or Polygons Only")
        break
del mxd
arcpy.AddMessage("\nDone!\n")


texts: Increase arrow symbol length by value in Arcmap

texts:
I use ArcMap to visualize surface displacements. I have a point shapefile with magnitude and angle values, which I plot as simple symbols with size=magnitude and rotation=angle.

As symbol I use arrows which length should increase at higher magnitudes. Anyone have an idea to realize that? My arrows increase in size but I didn't find a method for increasing the length only.

As you can see the size of arrows change with magnitude. I want that all arrows have the same size and only change in length.


texts: Create a new File Geodatabase(10.0) using data from existing Filegeodatabase(9.3)

texts:
I want top create a new file geo-database using python, which I am able to do using arcpy.CreateFileGDB_management(Fpath,"TestGDB","10.0")

But I want to transfer all data from an existing File geodatabse in ArcGIS 9.3 to this newly created geodatabase, what is the best way to do so using python?

I used  CopyFeatures_management(r"\Servername\Datafolder\Test.gdb","C:/Test/TEstGDB2.gdb")
but it gives me an error stating ERROR 000732: Input Features: Dataset\Servername\Datafolder\Test.gdb does not exist or is not supported


texts: Online GIS Masters / MBA / MS / MA with a focus on business

texts: Does anyone know of any GIS graduate programs that focus on the business elements of GIS? I am thinking of something like an MBA but with a GIS focus or a MA / MS with a business focus. Online only. Something like but GIS and Business Management. A web search is coming up with little.


texts: Looking for a fine U.S. Latitude and Longitude grid (0.25, 0.5 degree)

texts: There is a global latitude and longitude grid data available at 1 degree resolution and higher (http://www.arcgis.com/home/item.html?id=f11bcdc5d484400fa926dcce68de3df7). Is anyone aware of a finer resolution dataset (at 0.25 and 0.5 degree)? Or alternatively a way to make it in ArcGIS/GDAL/Python using the 1 degree product?


texts: Line first,last node distance along another line QGIS

texts:
I am using QGIS 2.18.6 . I have created two line layers. First layer "testroad" with long line (road). Second with short lines along,near road. I would like to know distance of first and last point of every line in second layer along line of first layer without creating another point layer. I've tried using expression in field calculator:

line_locate_point(geometry:=geometry(get_feature('testroad', 'id', '1')), point:=end_point($geometry))

But it always returns 0. I am new to QGIS and don't understand what is wrong with:

point:=end_point( $geometry)

rest of expression works.


texts:
Drawing a baseline in an attempt to make a grid. Have managed to draw baseline but cannot select using cad tools "select one linesegment" tool. Worked previously and I'm not aware that I'm doing anything different.

Baseline was made using this process:
Layer -> create layer -> new shapefile layer (line button checked, correct CRS) -> toggle editing -> add feature tool -> create two vector points to make a line -> Save Layer Edits

After this I try to select and it just doesn't register. Can select using the "select features by area or single click" tool but cannot create parallel lines using this (which is what I'm trying to do).

Possibly a bug, submitted a ticket for it at 13:05 GMT 09-07-15 will update if get response.


texts: Differentiating between stacked parcel shapefile and general parcel shapefile using Python

texts:
I am totally new to this domain. I would like to detect stacked parcel shapefile and mapped respective data or outlines created to all stacked parcels. In a general case where only one parcel on land it's easy to map one to one mapping. But when one building apartment has 20 floors and each floor has a different parcel for this case I want to map data to each parcel like one to many relationship.
Is it possible in Python to get some differentiable data between stacked parcel shapefile and general shapefile?
Stacked parcel example :


texts: QGIS merge multiple layers and keep highest valued polygon / grid cell

texts:
I'm still at a beginner level at QGIS. I'm aware of MMQGIS toolbox but I still can't figure out how to solve this problem.:
I got multiple layers indicating the number of events in a grid cells within a certain time. Figure 1 (orange) shows the number events in grid cells between 1AM and 4AM and figure 2 shows events between 4AM and 9AM.
Now I want to merge these layers so that the new layer will consist out of orange and blue grid cells depending on the higher number of events in the overlaying grid cell. If the number of events in a orange grid cell is higher than in a blue grid cell then the grid cell should be orange and vice versa.
This problem sounds simple quite simple to me though. Maybe you got some hints how to approach this problem.
Figure 1:Figure 2:

EDIT (my solution):
I did a quite long workaround with Python before where I merged all my data (which I saved in QGIS from .shp to .csv) with all timeframes (1AM-5AM, 5AM-9AM, 9AM- 1PM, 1PM-5PM etc.). I added a column with the appropriate timeframe. Then I grouped my data according to the field IDs. Finally, I calculated the maximum value within the group.
The output: I had the field ID with the highest value and the appropriate timeframe. I saved the data to CSV.
Afterwards, I merged my shapefiles with MMQGIS > Combine > Attributes Join from CSV file to add the timeframe attribute to the Field IDs. I saved to output layer and selected Categorized style with the timeframe column to set and display my values accordinly.


texts: How to use the replace string function with a path with special characters?

texts:
I am trying to replace a string path in a shape by replace function. I have this string:    C:\Users\me\Desktop\p1.jpg and I want to replace only thisC:\Users\Desktop\me\Desktop for 'C:\Users\main\peter and get C:\Users\main\peter\p1.jpg

I put this code replace('C:\Users\me\Desktop','C:\Users\me\Desktop','C:\Users\main\peter')but I getC:?sers?ain?eterwithoutp1.jpg` and with ? instead of \


texts:
I want to programmatically create a GeoJSON file using polygons from a shapefile but adding attributes from my own application.

This is easily done for a shapefile:

def create_data_dayer(self,varlist, data):
    """
    Creates a new shape to contain data about nodes.
    varlist is the list of fields names associated with
    the nodes.
    data is a list of lists whose first element is the geocode
    and the remaining elements are values of the fields, in the
    same order as they appear in varlist.
    """
    if os.path.exists(os.path.join(self.outdir,'Data.shp')):
        os.remove(os.path.join(self.outdir,'Data.shp'))
        os.remove(os.path.join(self.outdir,'Data.shx'))
        os.remove(os.path.join(self.outdir,'Data.dbf'))
    # Creates a new shape file to hold the data
    if not self.datasource:
        dsd = self.driver.CreateDataSource(os.path.join(self.outdir,'Data.shp'))
        self.datasource = dsd
        dl = dsd.CreateLayer("sim_results",geom_type=ogr.wkbPolygon)
    #Create the fields
    fi1 = ogr.FieldDefn("geocode",field_type=ogr.OFTInteger)
    dl.CreateField(fi1)
    for v in varlist:
        #print "creating data fields"
        fi = ogr.FieldDefn(v,field_type=ogr.OFTString)
        fi.SetPrecision(12)
        dl.CreateField(fi)

    #Add the features (points)
    for n,l in enumerate(data):
        #Iterate over the lines of the data matrix.
        gc = l[0]
        try:
            geom = self.geomdict[gc]
            if geom.GetGeometryType() != 3: continue
            #print geom.GetGeometryCount()
            fe = ogr.Feature(dl.GetLayerDefn())
            fe.SetField('geocode',gc)
            for v,d in zip (varlist,l[1:]):
                #print v,d
                fe.SetField(v,str(d))
            #Add the geometry
            #print "cloning geometry"
            clone = geom.Clone()
            #print geom
            #print "setting geometry"
            fe.SetGeometry(clone)
            #print "creating geom"
            dl.CreateFeature(fe)
        except: #Geocode not in polygon dictionary
            pass
        dl.SyncToDisk()

since I have all the geometries on a dictionary by geocode (self.geomdict) I simply create the features, set the fields and clone the geometries from pre-existing layer (code loading that layer omitted for simplicity). All I need now is a way to generate the GeoJSON from the combination of fields and geometries, naturally with the help of OGR to get the rest of the file right (CRS, etc. as from the source map)

How do export the feature collection generated as above?


texts: Can I replace a published vector layer with another one?

texts:
I have published a vector layer of polygons in GeoServer.
Can I replace it, with the same WMS and WFS services, with a newer layer?  It is basically the same layer with added features.


texts: Create a circle of defined radius around a point and then find the overlapping area on a shapefile

texts:
My use case is related to calculating the area of a polygon in a shapefile that is covered by another polygon.

Specifically, I want to do the following:

Create circles that have a 100m radius, from a data.frame of long/lat points
Overlay these circles onto a shapefile
Calculate the area of each of these polygons that are covered by these circles.
Ideally return back a dataframe where each row contains the area covered by all circles for each polygon.

Here I create some simulated data points.

mean_lat <- 46.07998
sd_lat <- 0.1609196
mean_long <- 8.931849
sd_long <-  0.1024659

dat_sim <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

Then we get the shapefiles

library(raster)
library(tidyverse)
library(sf)

# regular sp format
dat_ticino_sp <- getData(name = "GADM",
                         country = "CHE",
                         level = 3) %>%
  subset(.@data$NAME_1 == "Ticino")

# convert to sf format
dat_ticino_sf <- getData(name = "GADM",
                         country = "CHE",
                         level = 3) %>%
  # convert to simple features
  sf::st_as_sf() %>%
  # Filter down to Ticino
  dplyr::filter(NAME_1 == "Ticino")

So what I want to do is :

create circles with a 100m radius in the dat_sim, to overlay onto the shapefile dat_ticino_sf or dat_ticino_sp files.
For each polygon in the shapefile, calculate how much of that polygon is covered by circles.
Return a data frame where each row (polygon) contains the area covered by the circles.
Plot this.

Bonus points if this can be done in Simple Features.


texts: Using ArcPy: adding feature class and fields in an un-saved mxd

texts:
With ArcPy I have created a tool that creates an empty feature class with pre-defined fields and adds it to the TOC.  It works quite well, however if I use the tool to generate a feature class in an unsaved MXD it hangs.  I can get around it by naming and saving the MXD (it'S a bit cumbersome).
I'd like to change the code so that I can add a feature class to the TOC without necessarily creating an MXD project each time.

#------------------------------------------------------
# currently only works on a saved mxd.  needs to be able to add to an unsaved docuement.
#----------------------------------------------------------------
import arcpy
#import time feature - to insert date into field, but not functional yet.
from time import strftime
from arcpy import env
env.workspace = "CURRENT"

#--------------DEFINE VARIABLES---------------
in_wrkspace = arcpy.GetParameterAsText(0)
projectNumber = arcpy.GetParameterAsText(1)
FeatureDescription = arcpy.GetParameterAsText(2)
in_GeomType = arcpy.GetParameterAsText(3)
spatialRefCode = arcpy.GetParameterAsText(4)

#-----------------DEFINE EXPRESSIONS-----------------
FeatureClass = '_'+ str(projectNumber) + '_cm_' + str(FeatureDescription) + '_' + strftime("%y%m%d")

#Create feature class within workspace and add fields
FC = arcpy.CreateFeatureclass_management(in_wrkspace,FeatureClass,in_GeomType,"","DISABLED","DISABLED",spatialRefCode)
for FeatureClass in FC:

    # Process: classe
    arcpy.AddField_management(FeatureClass, "classe", "TEXT", "", "", "50", "", "NULLABLE", "NON_REQUIRED", "")

    # Process: objet
    arcpy.AddField_management(FeatureClass, "objet", "TEXT", "", "", "200", "", "NULLABLE", "NON_REQUIRED", "")
    #Process: source
    arcpy.AddField_management(FeatureClass, "source", "TEXT", "", "", "80", "", "NULLABLE", "NON_REQUIRED", "")
    #Process: date
    arcpy.AddField_management(FeatureClass, "date", "TEXT", "", "", "10", "", "NULLABLE", "NON_REQUIRED", "")
    #Process: remarque
    arcpy.AddField_management(FeatureClass, "remarque", "TEXT", "", "", "250", "", "NULLABLE", "NON_REQUIRED", "")
    #Process: traitment
    arcpy.AddField_management(FeatureClass, "traitment", "TEXT", "", "", "250", "", "NULLABLE", "NON_REQUIRED", "")
    #Process: no_seq
    arcpy.AddField_management(FeatureClass,"NO_SEQ", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
    #make feature layer

#Create a layer file to add to MXD
temp = '_'+ str(projectNumber) + '_cm_' + str(FeatureDescription) + '_' + strftime("%y%m%d")
arcpy.MakeFeatureLayer_management(FeatureClass,temp)

#-----------------add the layer to data frame--------------------------
#source: http://gis.stackexchange.com/questions/62896/how-to-add-feature-class-to-mxd-with-arcpy-python
#arcpy.MakeFeatureLayer_management("in_memory",FeatureClass)
MXD = arcpy.mapping.MapDocument("CURRENT")
DF = arcpy.mapping.ListDataFrames(MXD)[0]
Layer = arcpy.mapping.Layer(temp)
arcpy.mapping.AddLayer(DF,Layer,"AUTO_ARRANGE")
MXD.save()
#delete in memory cache layer

del MXD #, DF, Layer
arcpy.Delete_management(Layer)

Many thanks!


texts: Does Quantum GIS recognize QDGC information?

texts:
So a bit of background on what I am working on:
We are in the process of assisting a wildlife institution in South Africa with trying to take their information and surfacing it in a commercially available BI tool.

As part of their data they store Quarter Degree Grid Cell information about an extremely endangered mammal (Not Rhinos). They do not want to specify precisely where they spotted the animal due to the risk of poachers getting hold of the co-ordinates and poaching the animal.

The problem:
I have not come across a tool that naively recognizes QDGC information, and I don't want to start writing custom code for them as they don't have the skills to support the solution should something go wrong in future. This means that the solution that we will deliver will not be sustainable, which is definitely not ideal.

What I would like:
So my question is does Quantum GIS recognize QDGC information? If not, could the functionality be built in to the tool? If still not, could you please point me in the right direction?


texts: ArcMap: Combining County Polygons for 3 separate states into one layer with attributes

texts:
I have 3 separate layers with the counties for 3 separate states: WA, OR, & ID. I want to combine them into one layer "Counties" so they all share an attribute table with the same fields. I am having issues with the Spatial Join and Union tools where instead of joining the fields, it creates new fields with the same names. The following image shows what I am describing.

How can I either have the layers merge and their attributes combine into one field. OR could I create a new field and populate it using the field calculator? I've tried my luck at writing an expression but my SQL and Python are very weak.


texts:
Task: Display multiple symbologies (Gratuated symbols and Colours) on one layer

Problem: The symbology keeps overwriting itself, so I cannot get both symbologies at the same time

Methods Tried:

Export the graduated colour layer as .lyr and import. (Result:Overwrites)
Multiple attributes:(Result: Gives me all combinations, and that is not what I want)

I really don`t understand why I cannot get what I want to do in a few simple clicks?


texts:
I am trying to follow this Tutorial but in a different country. I am currently having a problem with cloud masks. I tried to visualize the different cloud masks of different patches at different dates. Sometimes the masks over-estimate the cloud coverage like the following screenshot, but in this case, I think it classified clouds and cloud shadow as clouds (not sure though):

But when it comes to small clouds they go undetected, which will most likely jeopardize later processing steps (crop type classification). Here is an example of that:

and this is the valid pixel count plot:

This is the time interval : time_interval = ['2019-01-01', '2019-12-31']
These cloud masks are generated using s2cloudless (a machine learning algorithm trained on S2-L1C images).
This is the code I used to download the images:
class SentinelHubValidData:
    """
    Combine Sen2Cor's classification map with `IS_DATA` to define a `VALID_DATA_SH` mask
    The SentinelHub's cloud mask is asumed to be found in eopatch.mask['CLM']
    """
    def __call__(self, eopatch):
        return np.logical_and(eopatch.mask['IS_DATA'].astype(np.bool),
                              np.logical_not(eopatch.mask['CLM'].astype(np.bool)))

class CountValid(EOTask):
    """
    The task counts number of valid observations in time-series and stores the results in the timeless mask.
    """
    def __init__(self, count_what, feature_name):
        self.what = count_what
        self.name = feature_name

    def execute(self, eopatch):
        eopatch.add_feature(FeatureType.MASK_TIMELESS, self.name, np.count_nonzero(eopatch.mask[self.what],axis=0))

        return eopatch

# Correspond to [coastal aerosol, B, G, R, vegetation red edge1,vegetation red edge2,vegetation red edge3, NIR, vegetation red edge4, water vapor, SWIR1, SWIR2] wavelengths BANDS-S2-L2A
band_names=['B01','B02','B03','B04','B05','B06','B07','B08','B8A','B09','B11','B12']

# TASK FOR BAND DATA and CLOUD MASK and CLOUD PROBABILITY
# add a request for S2 bands (downloading the data from the configurator)
# s2cloudless masks and probabilities are requested via additional data
add_data = SentinelHubInputTask(
    bands_feature=(FeatureType.DATA, 'BANDS'),
    bands = band_names,
    resolution=10,
    maxcc=0.1,
    time_difference=datetime.timedelta(minutes=120),
    data_source=DataSource.SENTINEL2_L2A,
    additional_data=[(FeatureType.MASK, 'dataMask', 'IS_DATA'),
                     (FeatureType.MASK, 'CLM'),
                     (FeatureType.DATA, 'CLP')])

# TASK FOR VALID MASK
# validate pixels using SentinelHub's cloud detection mask and region of acquisition
add_sh_valmask = AddValidDataMaskTask(SentinelHubValidData(),
                                      'VALID_DATA' # name of output mask
                                     )

# TASK FOR COUNTING VALID PIXELS
# count number of valid observations per pixel using valid data mask
count_val_sh = CountValid('VALID_DATA', # name of existing mask
                          'VALID_COUNT' # name of output scalar
                         )

I found out that fmask is also used for cloud masking but I couldn't download the python package in google colab.
So I want to know if there is any other way, using python, to mask those undetected small clouds.


texts: How to generalize a raster with noisy edges with ArcGIS Desktop?

texts: I'm trying to generalize a raster to get clean edges (see sample and "optimal" result below). However the borders bleed into each other to a rather high degree. I hacked away at the raster using majority filter, boundary clean and shrink, but it seems like I'm mostly shifting the noise around. What would be a good way to achieve a clean separation? The value range of the raster is 426-548, with 75 unique values. White = nodata.


texts: Show municipality in different colour than commune

texts:
I am using Atlas in Qgis 3.2 and would like to show municipality area in grey colour like in the picture (I made this using separate layer for a single municipality):

My goal is to be able to do the same but for the whole country with Atlas tool. I am already using an expression $id  = @atlas_featureid for showing only features contained in the Atlas current feature, which is a commune (not sure if this is the right word in English, it's a part of municipality).

I have a layer for administrative communes (which contains also names of municipalities) and layer for municipalities which I created from the first layer.

I think I should use an additional expression for only the municipality of the current feature to show in grey, but I can't figure it out.
Now all I can do is this (different municipality):


texts:
I'm trying to use the iteration function in ArcGIS model-builder. To avoid overwriting the output raster, I used %Name% inline variable. But the name of my files are getting bigger and bigger overtime (I use the outputs in another set of iterations). Is there a way to select only part of the file name and not all of it?

Thanks


texts: Will Grass GIS 6s atmospheric correction take raster as input for atmospheric parameters

texts:
I am working on atmospheric correction using 6s in Grass 6.4.4 version(i.atcorr module).

I want to use raster data set for atmospheric parameters. Presently I am creating an icnd (input condition)file having the atmospheric parameters in a format mentioned in grass help (please check https://grass.osgeo.org/grass64/manuals/i.atcorr.html).

But for accurate surface reflectance product one should use raster data (example like MODIS OMI-Aura_L3-OMTO3e_2014m0228_v003-2014m0304t144005.he5.nc )containing atmospheric parameters such as optical depth, water-vapor, ozone for any image.

Is it possible in Grass or do we have to use point data only for aerosols?


texts:
I changed the url i was using to fetch the esri js api to https://js.arcgis.com/3.7 and i can see in Fiddler that the call to that url is going out as ssl. But the file it pulls down (init.js) must be firing off other calls to fetch more components of the api and those are only going over http.
Short of hosting the api myself is there a way to fetch the whole enchilada over ssl?


texts: We are working on a QGIS layer contained in a project in .qgz format. So I did a deletion of the attribute table of a project layer. How can I recover the data that is not in the database backuped.


texts:
I am writing a tool within Python Toolbox that has two parameters:

param0 = arcpy.Parameter(
          name="a",
          displayName="A",
          datatype="GPFeatureLayer",
          parameterType="Required",
          direction="Input")

param1 = arcpy.Parameter(
          name="b",
          displayName="B",
          datatype="GPFeatureLayer",
          parameterType="Required",
          direction="Input")

params = [param0, param1]

I can´t figure out which if statement I should use in def updateMessages(self, parameters):. Here is a step-by-step what´s going wrong when I run the tool:

I choose a layer from the TOC as param0. My validation code is
if not parameters[0].hasBeenValidated:
    ... check e.g. if param0 has some attribute field ...
    ... it turns out that such field doesn´t exist, so...
    parameters[0].setErrorMessage("No field!")
The red error icon appears next to its label with the message "No field!" - everything is working as expected.
But then when I choose the second layer from the TOC as param1 (it has no validation code) the param0´s red icon disappears like there is no error with the param0! I am sure that the problem is with if not parameters[0].hasBeenValidated: as I don´t quite understand how it works. I tested several code variants with .hasBeenValidated, .altered etc. but with no success.

To be short:

I would need the red error icon to be displayed as long as there is an error with the param0 no matter what I do with other parameters and
(if possible - not mentioned above) ensure that param0´s validation performs only when the user changes the value of this parameter.

Could you give me a hint how to accomplish this?


texts: How to reference foreign keys in OGC API Features (pygeoapi)?

texts: I'm using a modified version of the PostgreSQL provider in pygeoapi, and I want to keep the relations between tables in my OGC API Features. Does the standard have a solution for this? I thought the solution was to add the reference table to the links, but I don't see a "relation name" in the list for this case.


texts: QGIS on Android: Access files on Google Drive

texts:
I have a set of shape files that I keep updated on Google Drive.  I have just installed QGIS on an Android tablet but I can not figure out how to get access to the google drive from the file menu.

Any ideas?

Russell


texts:
Using QGIS 3, I'd like to split a polygon using a line. My inputs are .geojson files, but that should be irrelevant:

Line:

{
"type": "FeatureCollection",
"name": "line",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "LineString", "coordinates": [ [ 3.54692, 49.76417 ], [ 3.88975, 49.22462 ], [ 4.64133, 49.42229 ], [ 4.58858, 49.06073 ], [ 5.18193, 49.21601 ], [ 5.4852, 49.06937 ], [ 5.41927, 49.33644 ], [ 5.11601, 49.43086 ], [ 5.31379, 49.61916 ], [ 5.88077, 49.44801 ], [ 5.99944, 49.27626 ], [ 5.78847, 49.02616 ], [ 5.36653, 48.86162 ], [ 4.78637, 48.68782 ], [ 5.07645, 48.49596 ], [ 6.60598, 48.29459 ] ] } }
]
}

Polygon:

{
"type": "FeatureCollection",
"name": "poly",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 4.49628, 49.90025 ], [ 2.66349, 48.93096 ], [ 4.93141, 48.08362 ], [ 6.43456, 48.69653 ], [ 6.34227, 49.92572 ], [ 4.49628, 49.90025 ] ] ] } }
]
}

Here is my code:

poly = QgsVectorLayer(r"D:\SplitPoly\poly.geojson", "poly", "ogr")
line = QgsVectorLayer(r"D:\SplitPoly\poly.geojson", "line", "ogr")

polyGeom = poly.getFeature(0).geometry()
lineGeom = line.getFeature(0).geometry()

splitGeom = polyGeom.splitGeometry(line.getFeature(0), True)

print('Result:')
print(splitGeom)

You can also swap the 'splitGeom line' with:

splitGeom = polyGeom.splitGeometry(lineGeom.asPolyline(), True)

But either way the result is:

(1002, [], [])

Any suggestions what is wrong with my script?


texts: Sampling tiles in regular grid with Earth Engine

texts:
I need to extract tiles of size 256x256 from an ee.Image over regions where I have labeled data. I understand how to take a random sample within polygons from an object returned from ee.Image.neighborhoodToArray() by using sample(), but I need a uniform sample over a grid defined by the extent of the labels. The image below shows an example of a label raster overlaid with tiles that I need to extract.

I understand I could do it by generating points in a uniform grid over my label shapefile and sampling at those points, but this seems cumbersome. Is there functionality in Earth Engine to do this operation?


texts: Vector differencing fails in QGIS (due to geometric errors)

texts:
I was trying to correct the country shapefile of Uzbekistan for Gas flares:

Input layer (the UZB_adm1) differenced by this layer.

In QGIS the result was that the overlapping regions - where the differencing should have taken place - were gone. I presumed it was due to wrong referencing of the UZB shapefile. Couldn't resolve it somehow.
After several desperate hours I went to gvSIG 2.2, did the differencing job there. And it worked like a charme.
Has anyone an Idea what was wrong in QGIS?


texts:
I am working with qgis api c++ 2.8.2, in c++/QT 4.5
I try to save my List of QgsField but I have the error message :
no matching for call to 'QgsVectorFileWriter'
I took the sample from here :

http://hub.qgis.org/issues/7223

Here is my code in C++

QList <QgsField> lMyAttribField;
QgsVectorFileWriter * writer;
writer = QgsVectorFileWriter("/tmp/myshapes.shp", "CP1250", lMyAttribField, QGis::WKBPoint, "0", "ESRI Shapefile");

I have  the same eror if I delete the * in the declaration of writer

I don't understand what's i forget ?


texts: GetLegendGraphic request string is too long?

texts:
I am setting the SLD Body of a WMS layer on the client side. I set the style in JSON and I have some functions that convert the JSON to SLD (OpenLayers functions). It works fine for displaying the layer with all its colors.

But when I use a legend (GetLegendGraphic from GeoServer), I have a limit in the size of the SLD Body. The limit of characters is about 20.000. When I perform the GetLegendGraphic request and the url is longer than (about) 20.000 characters, an error occurs and the legend is not displayed. BUT the layer is still correctly displayed (as you can see in the images below).

Below the 20.000 characters, the legend works fine.

I find it strange that the layer is correctly displayed and the legend isn't. Both are fetched via GeoServer requests.

Could you help me solve this?

I wanted a way to accept long strings in the GetLegendGraphic or a way to compress the string with the SLD body (so that it has fewer characters).

EXAMPLES:

SLD Body with 18006 characters (this shows a correct legend)

SLD Body with 24245 characters (this does not show the legend, still draws the layer correctly)


texts: How to combine several datapoints in one location to a single point varying in size accordingly?

texts:
This might be a very simple thing, but I just can't figure it out how to do it in QGIS (alternatively I'm also grateful for explanations for ArcGIS).
I have a fisheries dataset with multiple catches in one location. I'd like to display these catches as a dot varying in size according to the number of catches - like dots displaying population size in a city. However, the excel sheet I have is one catch each row, with the same coordinates but different attributes (length, weight, etc). Is there a tool in QGIS that can merge all points in one location (other attributes not needed then)?
Thanks a lot for your help!


texts: Most timely source of California Fire Map data?

texts: Fire season always has me scouring the web for satellite fire map information. Unfortunately many of the locations I find are outdated or not maintained well. What would be the best site to go to for the most timely fire data?


texts: Additional features added when using "join attributes by location" QGIS

texts:
I'm getting some weird output after performing the "Join attributes by location" tool inside QGIS 3.10. Below is a shapefile representing watersheds in South Africa. I'm only interested in the watersheds in direct contact with the sea (in green). As you can see, displayed in white are watersheds I rejected as they are endorheic basins, i.e. with no outflow.

I also have a point shapefile, each point is an average rainfall.

What I'm trying to do is very simple, it is just to join the rainfall data to the watershed polygons, so I can have an idea of the rainfall amount in each one. However, when I perform the tool, rainfall data is assigned to every basin, even the endorheic ones, as shown in the picture below.

The join tool was simultaneously performed using either "intersects", "contains" and "within" as geometric predicate. My goal is just to have the rainfall data for my targeted basins, not the endorheic ones.


texts: Geoprocessing all shapefiles of folder using ArcGIS ModelBuilder?

texts:
I have 200 shapefile in a folder. All of them are similar. I need to select by attribute then apply the Eliminate tool and save the shapefiles.

How do I do this automatically with all the feature classes using ModelBuilder?