texts: Converting NetCDF to Grib2 in Python on native Windows. Which package?

texts:
Problem statement:
I need to implement an easy-to-share script compatible with Windows environment that allows me to open a NetCDF, manipulate it, and save it in Grib2 format.
Attempt 1: iris / iris-grib
I initially implemented my code with success using iris and iris-grib in a Linux environment. However, iris-grib depends on cf-units/udunits2/grib_api packages, which have serious compatibility issues with Windows (e.g. [1]), so it would be hard to make the code easy-to-share.
Attempt 2: xarray + cfgrib
Later, I tried to perform something like:
import xarray as xr
import cfgrib

my_ds = xr.open_dataset(INPUT_FILE_PATH)

[...some manipulation on my_ds...]

cfgrib.to_grib(my_ds, OUTPUT_FILE_PATH, grib_keys={'edition': 2})

Which worked great until the cfgrib.to_grib(...) line. There I got the error:
ValueError: required grib_key 'gridType' not passed nor auto-detected.
As my input NetCDF file is in EPSG:4326 projection (WGS-84), I tried to set the grib_keys arguments as grib_keys={'edition': 2, 'gridType': 'regular_ll'}.
ValueError: sample not found: 'regular_ll_sfc_grib2'
Which seems to be an issue with the cfgrib installation (done using conda) not being able to find its own ancillary internal files.
Attempt 3: cdo
Ok, so I saw that there is a Python package that wrappers cdo. Great! So with some file writings/readings, I can do what I want. Right?
Not really.
The cdo package installed with conda install cdo is indeed just a wrapper, so that the system is expected to have the own cdo command available and ready to run. As running cdo on Windows requires a cygwin set up, making the code easy-to-share would be somewhat really messy.
So...
What is a good and easy-to-share way to save a grib2 file on Windows 10 out of a NetCDF file?
[1] https://stackoverflow.com/questions/38979819/iris-for-anaconda-fails-to-read-grib


texts: Inserting landsat image WMS to MapServer map file?

texts:
I am trying to insert this landsat image to my Mapserver map file. However, I can not manage it. Does someone have an idea?

Landsat image WMS information:

Landsat from onearth.jpl.nasa.gov
Bounding box: (-180.0, -90.0, 180.0, 90.0) (show on OSM)
Projection: NASA:4326
WMS half-link: http://irs.gis-lab.info/?layers=landsat&
Tiles URL: http://irs.gis-lab.info/?layers=landsat&request=GetTile&z=!&x=!&y=!

my mapfile layer:

LAYER
    NAME global_landsat
    TYPE RASTER
    STATUS DEFAULT
    CONNECTIONTYPE WMS
    CONNECTION "http://irs.gis-lab.info/?layers=landsat&"
    MINSCALE 20000
        METADATA
            "wms_server_vesion" "1.1.1"
            "wms_srs"  "EPSG:4326"
            "wms_format" "image/png"
            "wms_name" "modis,global_mosaic"
        END
    END


texts: Qgis 1.9 on Ubuntu 12.04 constant crashes as soon as I try to do anything

texts: I can get a google earth physical layer open but as soon as I try to open a vector file it crashes-just goes away-no error message, nothing


texts: Why does the attribute Table not maintain input data?

texts: for a longer time I have a problem with the attribute table. I create the polygons on the map, I fill up the table with attributes. Then, if I click "save", all the data from the table are being deleted and in all fields of the table a word "error" does appear. After the restart of the program, only that columns containing the number values are being kept. That ones which contain the words are deleted and the word "null" does appear in all the fields. Is somebody able to help me?


texts: Sum every points value layer that intersect line and automatically update label sum if more points intersected

texts:
I am using 3.12 QGIS version. I have got two point layers, and one line layer.

I wanna put a label on line layer, that calculate "sum" of intersected points or any mathematical calculation.

The main purpose is collecting any points (value) which intersected with line, but to keep in mind, we assume there are thousands point intersected, so when we build the code in expression dialog, we dont need to define the points layer, because generalizing the solution would be better.


texts: How to change the format output from ftl

texts:
I have have a view in GeoServer from PostgreSQL and I like the default free marker template output in Layer preview but I would like to change it to where to have column names of Field name and Field Value.

Since the default content.ftl is iterating through and outputting all the content horizontally is there a easy way to make the table print vertically?


texts: Applying Edge Correction in ArcGIS Average Nearest Neighbor?

texts:
I am using ArcGIS 10.3 with the tool Average Nearest Neighbor from the Spatial Statistics toolbox to describe the data. I want to avoid mistakes caused by Edge Effect.

In the Ripleys K Function tool there is an option to control that (there are a few options of what to do with Edge Effect correction) but not in Average Nearest Neighbor.

How can I avoid Edge Effect mistakes while using Average Nearest Neighbor tool?


texts: Distributing processed SRTM data legally?

texts:
I'm working on different maps, based on SRTM data.

The goal is to lay out the maps at a park entrance against a little money for the printing and my personal costs. The SRTM 4 data from CGIAR I origanally used prohibits the commercial use, so I'm not sure, if I'm allowed to do so.

Is there any SRTM data, free to use for that kind of purpose?


texts: Conditionally adding markers to map in R leaflet and shiny

texts:
I am trying to add some markers to a leaflet map in R based on conditions that I select from an selectInput command. The selectInput command lets the user choose three different numbers (1, 2 or 4). Based on this number the respective marker should be added to the leaflet map. In the end I am attempting to only add one marker based on the selected number of selectInput.

My problem: I am trying to use ifelse() statements, but I get an error that the argument could not be interpreted as a logical condition. I am assuming that this might be linked to the piping of my leaflet command? My code with some sample data is below.

library(shiny)

# create some sample data containing lat/lon coordinates
data1 <- as.data.frame(cbind(c(46.5728, 46.5729, 46.5730, 46.5731),c(6.9385,6.9386,6.9387, 6.9388)))
colnames(data1) <- c("Lat", "Lon")
data2 <- as.data.frame(cbind(c(46.4728, 46.4729, 46.4730, 46.4731),c(6.8385,6.8386,6.8387, 6.8388)))
colnames(data2) <- c("Lat", "Lon")
data4 <- as.data.frame(cbind(c(46.3728, 46.3729, 46.3730, 46.3731),c(6.7385,6.7386,6.7387, 6.7388)))
colnames(data2) <- c("Lat", "Lon")

ui <- fluidPage(
  tabsetPanel(tabPanel("TEST PANEL",
                       fluidRow(leafletOutput("Map", height = 500),
                       wellPanel(fluidRow(selectInput("Number", "Select number:", c(1,2,4))))))))

server <- function(input,output, session){
  output$Map <- renderLeaflet({
    m <- leaflet() %>%
      addProviderTiles('Esri.WorldImagery', options = providerTileOptions(maxNativeZoom=19,maxZoom=100)) %>%
      if(input$map == "1"){
        addCircleMarkers(data = data1,
                       lat =  data1$Lat,
                       lng =  data1$Lon,
                       radius = 4, col = rgb(215,48,39, max = 255))
      }
      if(input$map == "2"){
      addCircleMarkers(data = data2,
                       lat =  data2$Lat,
                       lng =  data2$Lon,
                       radius = 4, col = rgb(233,25,78, max = 255))
      }
      if(input$map == "4"){
      addCircleMarkers(data = data4,
                       lat =  data4$Lat,
                       lng =  data4$Lon,
                       radius = 4, col = rgb(23,25,78, max = 255))
      }
      addScaleBar(position = "bottomleft", options = scaleBarOptions()) %>%
      leaflet::addLegend("topright",
                         colors =c(rgb(215,48,39, max = 255),  rgb(233,25,78, max = 255), rgb(23,25,78, max = 255)),
                         labels= c("No 1", "No 2","No 4"),
                         title= "Number",
                         opacity = 1) %>%
      addMeasure(primaryLengthUnit = "meters", decPoint = ".")
    m
  })
}

shinyApp(ui, server)


texts: Using iface.showAttributeTable() on QgsVectorLayer retrieved by getObject()

texts:
I'm unsure if i'm properly using the processing.getObject() function. I'm creating a python processing script for QGIS 2.0 and 2.2, but can't seem to get an attribute table to open programmatically. The code below works fine in the Python Console, but causes QGIS to crash when I tried as a Processing Script. The only difference was that the layer is called by QgsVectorLayer() in the Python Console and processing.getObject() in the Processing script, thus my reasoning behind my question.

from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from processing.core.VectorWriter import VectorWriter
import numpy
import processing
import csv
import qgis.utils

#Here we define the input and outputs
#====================================
##[Example scripts]=group
##Couche_compteurs_individuels=vector
##Couche_compteurs_generaux=vector
##Couche_canalisations_sur_compteur=vector

#And here is the body of the algorithm
#=======================================

#input layers values are always a string with its location.
#That string can be converted into a QGIS object (a QgsVectorLayer in this case))
#using the processing.getobject() method
indivLayer = processing.getobject(Couche_compteurs_individuels)
grosLayer = processing.getobject(Couche_compteurs_generaux)
resLayer = processing.getobject(Couche_canalisations_sur_compteur)

CG_IDlist = []

def AEPstats():
    Gros_Feat = grosLayer.getFeatures()
    grosIDselection = []
    grosidNull = False
    for feat in Gros_Feat:
        ID_idx = grosLayer.fieldNameIndex('ID_CG')
        #selection list of null values
        if feat.attributes()[ID_idx] == NULL:
            grosIDselection.append(feat.id())
            grosLayer.setSelectedFeatures(grosIDselection)
            grosidNull = True
        #Putting existing attributes in list
        else:
            CG_IDlist.append(feat.attributes()[ID_idx])
    #Open msg box, attribute table with selected items, and end function
    if grosidNull == True:
        QMessageBox.warning(
            None,
            'Probleme de saisie',
            "Completez identifiant(s) de reference ('CG_ID') de Compteurs Individuel")
        qgis.utils.iface.showAttributeTable(grosLayer)
        return

Thanks in advance


texts: How to adjust the Record Dates field in ArcGIS 10.3.1 for Desktop?

texts:
When I use the Editor Tracking in polygon layer, the Date time that been created is not true:

and my real time was:

I read:

http://desktop.arcgis.com/en/arcmap/10.3/manage-data/editing-fundamentals/enabling-editor-tracking.htm
In arcmap how to modify time value in editor tracking fields

but didn't understand how to change the Record Dates field in the Editor Tracking tool in the data management tools-any time,after I add the layer to the tool, the Record Dates option is blocked to only the UTC option:

I also worked on another polygon layer that don't have ant fields of Editor Tracking, and it still doesn't worked.


texts:
Currently there is the following software in use:
Win10 Clients, MS Access, ArcView 3.2 from year 1999

There are 3 common workflows, I describe one:

Within the MS Access Form a Dataset is selected
Then I click the „Show me in GIS“ Button in MS Access Form.
MS Access opens ArcView and an ArcView Project with a DDE Connection, not only is this
DDE Connection opening ArcView, it also takes the ID of the selected Dataset and runs a script within ArcView.
The script in ArcView gets the ID from the DDE connection, goes to a known shapefile and looks into the attibute table if the ID is present, if yes, it zooms to the proper Geometry.

So no Access Database is directly connected to ArcView, but only some data is sent and received.

The question:
How it is possible to achive the same workflow with QGIS? (The MS Access Database will stay literally forever.) I can't use a DDE connection with QGIS as far as I know.

I am not asking here for a perfect solution. A starting point would be good (PyQGIS, batch-file, ...)


texts:
i saw your code to create a buffer around a point.
on the link:
How to add a buffer on a point in OpenLayers 3

in following code..

var pointFeature = new ol.Feature(
new ol.geom.Point([0,0])
);

How can I pass the coordinate of a point where a user click on map.
I have 2 layer in my map which holds the different features with in it.
I want to create buffer on a particular layer which will be given by user.
I have one popUp implemented on map.onclick already.
And somehow i have implemented buffer also on map's click event(separately).

How can I manage both popUp and Buffer on map's click event?


texts: I drew the wind barbs on the map based on the coordinates data with the step of 2 degrees. But I am not satisfied with the fact that after the map reprojecting the size of the wind barbs changes. Is it possible to draw a wind barbs with exactly the same size with GeoTools?


texts: Problem getting correct area for polygon and choosing a CRS

texts:
Using GeoPandas and I'd like to calculate the areas for the geometries in my dataframe.  I'm also doing other things like finding the area of overlapping regions of geometries from different dataframes, but I think an answer to the basic question will solve all my problems in area calculations.
The polygon coordinates are all saved as "naive geometries" and I've set the CRS to 4326.
To calculate areas, I only need a CRS that works reasonably well across the Tokyo area.
I have data from the government that includes the land area of administrative units, and using CRS=3095 gives very close results (so does {'proj':'cea'}).
One of my datasets is a hexagonal grid, where the side radius (apothem) is 125 meters, so the area of each hexagon is 54,127 m2.  However, when I calculate the area of the hexagon geometries using CRS=3095 I get 3,764,267,156 (which should be meters sq because the CRS documentation says it is in meters).
Using CRS=3857 I get 5,764,828,070.  That's very different, but still in the same ballpark as the other CRS compared to the known area of 54,127 m2.  Those are insanely large (and wrong) numbers for my hexes.
So my question is not so much choosing a CRS for my region, but rather how do I interpret the values coming from these calculations?  Knowing the area is 54,127 m2, how can set a CRS and perform the calculation of area on my geometries such that it will produce that value (to a reasonable approximation)?


texts: How to transform SRID 3035 to SRID 4326?

texts:
I have data in the ETRS 89 LAEA format (SRID: 3035). But I have problems to transform it.

Example data:
    N2736 E4541

I want to transform it to WGS84 (SRID: 4326), to check the coordinates on a map.

My transformation looks like this:
    SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(27.36 45.41)',3035),4326))

I get: POINT(-29.0867265338903 12.9940786314927)
But the point should be in Germany and not in the atlantic ocean :(


texts: How to highlight a clicked Feature Layer using ArcGIS JS API 4.x

texts:
Esri recently released version 4.0 of the ArcGIS Server JavaScript API, which represents a substantial re-write compared to version 3.x

Using version 3.x, a feature layer supported a long list of events, including click and mouse-over. For example, it was easy to highlight the boundary of a polygon by listening to the feature layer's mouse-over event and accessing its geometry, as in this sample.

Using version 4.0, the functionality matrix says that Mouse Events are "not planned" and that we should use View Events instead.

How can you highlight the clicked feature on a feature layer, from a view's click or hover?


texts: Grid (Moving average) how to set moving distance?

texts:
I'm using Grid (Moving average) to get mean values of Point attribute in a circle of 20 m radius.

What is determining the step distance / the output pixel size? I can't find it in the documentation. I would like to set I manually to for example 10 m.


texts:
I somehow closed the layer TOC in QGIS 2.0.1.
How do I get it back? Feeling a bit lost here.


texts: Polygon boundary is not properly connected in R

texts:
I have two points sets, when plotted using geom_point, they look like:

What I want to do is draw a polygon using these two datasets, and I expect to get a ring-shape polygon. However, when I use geom_path to draw the boundary, I got a weird result that boundaries are not properly closed.

I think the problem may due to the order of the dataframes, but I cannot figure it out. Any suggestion?


texts:
I would like to write a Python expression function which calculates the average of columns (variable length parameters), skipping NULL values.

@qgsfunction(args=-1, group='Custom')
def col_avg(vals, feature, parent):
    v = [ x for x in vals if x is not NULL]
    if len(v):
        return sum(v) / len(v)
    return 0

If there is a NULL value among the input parameters the result is NULL, too.
If I use an expression in Select by Expression for example:

my_col is not NULL

it works.
I use QGIS 2.18.7 on Ubuntu 16.04

* addition *
Here is an image, how I call my function:


texts: How can I find out which  postgis.backend (geos or sfcgal) I have?


texts: Embedding feature chart/graph in pop-up window in CartoDB?

texts: I would like to be able to embed a chart/graph of attributes for each feature within a webmap. I've had a look around and haven't found anything definitive on the turorial pages of how to achieve this.


texts:
I am developing a car-sharing application and I want to check if the approach I have is correct. I am using Postgres with PostGIS.

I have a list of routes that users publish, where they indicate that they typically commute from a point A to a point B (i.e. I am going from Milan to Venezia). All those routes are saved in the database as polylines that represent user journeys.

Then I need to search all those routes giving an input route (the red line in the picture), looking for routes that intersect my input or are close by (let's say 10 km deviation).

My solution is to convert everything to linestrings and use the Hausdorff distance to get close paths. Since I am fairly new to this approach, before proceeding further I have some doubt:

Is this approach computationally feasible?
Does it take into account the direction of the two polylines?
Is it a good approximation of "similar" routes?
Are there optimizations I can do to improve the overall search?


texts:
Which way do I calculate bottomlines of valleys or even more preceise:
which way do I "find" valleys (for example suggested tectonic faults).

Sample below shows a height profile (profile from line) based on Aster
data and it's that I at most want to follow the larger valleys (long
arrowed lines).

I could go for hydrological methods (even though havent been working
with them), but using hydro plugins/tools a small breakthrough to lower areas would cause a runoff, while I want to follow the big valley.

Any ideas?

--- (edited)
To underline my problem, I've just taken a circa bottom line profile
of the valley pointed to by the most left arrow in pic-1 (arrow'ed).
The intersection is at about 8km of 1st diagramm (with arrows) and at
km 35 of second profile.

---- (edited)
Just to spread some ideas:
Rotating the 1st image arount the center between point 1 an 2 (1st and 2nd arrow, Height Profile) also shows "valleys" with different idth/height ratios, but still significant.

QGIS 1.8, Ubuntu, Aster Elev Data

First and second "sharp" valley from left (large arrows) are of most interesting to get the idea.
It might help to focus on focusing on strong beeing above/below average. Just an idea :-/


texts: is it possible to show all labels of WMS layer at any extent without SLD ol5?

texts:
I am showing WMS layer that shows label on it by defining SLD in GeoServer. The problem is that is is not showing all labels at particular extent.

But when I zoom in, all labels are displayed.

My problem is that I am exporting the map at some extent then all label is not showing when I see the PDF.

Is it possible to show all labels in map PDF?


texts: setting a PictureFillSymbol via ArcObjects not wanted/allowed?

texts:
Have a problem in setting a PictureFillSymbol in polygon-layer-properties via ArcObjects. I just create the Bitmap not from file, it is "painted" in runtime. After loading the layer and setting the properties, the symbol in legend and map is correct shown, even colors. But in layer-properties-dialog the foregroundcolor could not be set, it is disabled/greyed. Also the checkbox for swap the colors. It is a 1bppindexed BMP, just don't know, why ArcGIS this could not see.

Here is a snippet from my code setting the symbol (Maybe somone has an idea?):

System.Drawing.Image img = CreateBitmapFromStipple(stipple);
// Farbpalette setzen
System.Drawing.Imaging.ColorPalette cplr = img.Palette;
System.Drawing.Color colFG=System.Drawing.Color.FromArgb(255,col.Red,col.Green,col.Blue);
System.Drawing.Color colBG = System.Drawing.Color.FromArgb(255, bgcol.Red, bgcol.Green, bgcol.Blue);
cplr.Entries.SetValue(colFG, 0);
cplr.Entries.SetValue(colBG, 1);
img.Palette = cplr;

ISymbol fSymbol = new PictureFillSymbolClass();
(fSymbol as IPictureFillSymbol).Color = col;
(fSymbol as IPictureFillSymbol).BackgroundColor = bgcol;

stdole.IPictureDisp pd=VB6CompatibilityConverter.ImageToPictureDisp(img);
(fSymbol as IPictureFillSymbol).Picture = pd;

//(fSymbol as IPictureFillSymbol).BitmapTransparencyColor = bgcol;
(fSymbol as IPictureFillSymbol).SwapForeGroundBackGroundColor = textReader.GetAttribute("SwapForeGroundBackGroundColor").ToLower() == "true";
if (double.TryParse(textReader.GetAttribute("YScale").Replace('.', ','), out dHlp)) (fSymbol as IPictureFillSymbol).YScale = dHlp;
if (double.TryParse(textReader.GetAttribute("XScale").Replace('.', ','), out dHlp)) (fSymbol as IPictureFillSymbol).XScale = dHlp;
if (sOutLine != null) fSymbol.Outline = sOutLine;

  Added:

Maybe it is not allowed to set the Picture directly through Arcobjects?
In API reference I found out, picture is only a read-member, not write.

But however the set-Property is defined.


texts:
I made two river bank lines and I want to associate them.

I want to make a centerline using two rows as a base. I managed to make the centre line with a lot of work creating several lines from one margin to the other of the river and then using the centroid, but this solution is very tiring and I have to do this in several areas.


texts: Any way to remove 2 point polylines in bulk?

texts:
QGIS:
I have successfully created contours from a DEM and subsequently cleaned them up using generalizer tools in the following order:
1.  Vertex reduction
2.  McMinimans distance algorithm.

The contours look great, but I am stuck with millions of dots that seem to be 2 point polylines with both vertex's in the same location, making them appear as single points.

These points are not selectable using the selection tool either which is maddening.

Can anyone help me with an idea of how to delete them in bulk or perhaps, how to not get them in the first place?


texts: pyqgis: select features inside a buffer, then return some attributes of the selection as a list

texts:
I'm a python and pyqgis beginner (2 months-self taught) and I found it really hard to learn pyqgis API in comparison to the other APIs like kivy or Pyqt4.

I'm writing a qgis plugin that

user select 2 input layers and fields via combo box (which is fine).
creates a buffer around the first layer's features.(it's fine too)
select elements of the second layer that are inside the buffer, which user will choose the second layer via combo box also(here is the problem)
it has a couple more combo boxes that show the second layer's fields. and I want to save user-defined attributes of the second layer's selected features into a list of list.

here is the code:

def run(self):
"""Run method """
    # region adding vector layers
    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        if layer.type() == qgis.core.QgsMapLayer.VectorLayer:
            layer_list.append(layer.name())
        else:
            continue
    # endregion

    # region first layer inputs
    self.cInput = self.dlg.comboBox_input_Crime
    self.cInput.clear()
    self.cInput.addItems(layer_list)
    # endregion

    # region second layer, layer and combo boxes
    self.pgInput = self.dlg.comboBox_input_polygon
    self.pgInput.clear()
    self.pgInput.addItems(layer_list)

    # these are the target fields of selected features that must be written in a list of list
    self.pgFactor = self.dlg.comboBox__factor_polygon
    self.pgLyr = self.dlg.comboBox_layer_polygon

    # update fields
    def pgField_select():
        self.pgFactor.clear()
        self.pgLyr.clear()
        selectedLayerIndex = self.pgInput.currentIndex()
        selectedLayer = layers[selectedLayerIndex]
        fields = [field.name() for field in selectedLayer.pendingFields()]

        self.pgFactor.addItems(fields)
        self.pgLyr.addItems(fields)

    self.pgInput.currentIndexChanged.connect(pgField_select)
    # endregion

    # show the dialog
    self.dlg.show()

    result = self.dlg.exec_()

    if result:

        # get points as list
        self.pic = self.cInput.currentIndex()

        self.buff =  processing.runalg('qgis:fixeddistancebuffer', layers[self.pic], 10 , 5, True, None)

        self.buffLyr = qgis.core.QgsVectorLayer(self.buff['OUTPUT'], "buffer1", "ogr")

        qgis.core.QgsMapLayerRegistry.instance().addMapLayer(self.buffLyr)

        #here is the problem, it does not selects anything
        processing.runalg("qgis:selectbylocation", layers[self.pgInput.currentIndex()], self.buffLyr, ['touches'], 0)

however, I faced a nightmare learning pyqgis and completing this task.

how can I select features by location and save their attributes in a list of list like: [['river','Dallas'],['pond', 'Houston'],...]


texts: Export Magnetfield collected data to shapefile

texts:
I am currently using a Sokkia GHX2 to collect field data with Magnetfield. We can easily exchange collected data to a dwg file and there is an option to exchange to shapefile. However, when I try to export to the file to an ESRI - it creates an individual shapefile for each point.
Is there a way I can export the entire job and have the codes and names connected to each point?


texts: WorldWind wms services show reversed coordinate in ArcGIS Desktop

texts: I'm adding WorldWind wms services (http://data.worldwind.arc.nasa.gov/wms?) to ArcMap (ArcGIS 10.2.2) but maps are misplaced, looks like they have reversed coordinates (as explained here: http://viswaug.wordpress.com/2009/03/15/reversed-co-ordinate-axis-order-for-epsg4326-vs-crs84-when-requesting-wms-130-images/). Is there any workaround to get maps properly placed?


texts: Why can't I export XY data to a shapefile?

texts:
I have added a table from Access to ArcMap, and then used 'Display XY data' to plot points. When I then try to export the Events to a shapefile I get an error message telling me:

"Error exporting data.
Create output feature class failed"

Followed by:

"Warning: Some records may not have been exported correctly."

and then:

"There was an error exporting the layer."

Can anyone tell me why this is happening, and how I can export my data?

I was able to export an almost identical file this morning with no problems. I'm saving straight to my C drive, so no unusual characters in the files path. All the fields are either text or Long types, and none have long names/unusual characters in their names.

I've tried restarting ArcMap and my machine, and exporting to several different locations.

I'm using version 10, with Access 2007, and there are 7332 points in the table.


texts: Need to set the display expression for multiple layers in ArcMap

texts:
I have multiple layers in ArcMap that I split (using ModelBuilder, feature selector with a feature class to feature class) from a single layer, all within the same fileGDB.
I can use the batch "Apply Symbology from Layer" tool to match the symbology from the original layer. But I can't figure out how to set the correct field for the display expression, and it takes too long to go through and set each layer manually (~20+ layers).
Any thoughts?


texts: How to get polygons contained within another polygon, including those that touch the border in BigQuery Standard SQL?

texts:
I'm trying to get blockgroups contains within a designated market area polygon in BigQuery. I've tried using st_contains and st_covers but I still only get the completely contained ones (and not the bordering ones:
SELECT a.blockgroup_geom as the_geom,
a.geo_id,
c.total_pop

FROM `bigquery-public-data.geo_census_blockgroups.us_blockgroups_national` a

join `bigquery-public-data.census_bureau_acs.blockgroup_2018_5yr` c

on a.geo_id = c.geo_id

join `bigquery-public-data.geo_us_boundaries.designated_market_area` b
on st_covers(b.dma_geom,a.blockgroup_geom)
where b.dma_name = 'Milwaukee, WI'

What do I need to do to get all of the blockgroups within the dma polygon?


texts: Convex hull area: what are the output units?

texts:
I have what I hope is a simple question. I'm doing convex hull calculations. What I have is area(convex_hull($geometry)).

I get a number and that's all well and fine. Since I need to get a ratio of that area with the original polygon, I need to know what units the convex hull area is in. Does anyone have an idea on this?

Also, am I doing something naive?


texts: PySAL, ValueError: cannot assign slice from input of different size

texts:
I want to calculate local moran's I statistics using pysal with Queen contiguity but I'm getting error.
import geopandas as gpd
from pysal.explore import esda
from pysal.lib import weights
gdf = gpd.read_file("\\CA_Places\\CA_Places_TIGER2016.shp")
w_queen = weights.Queen.from_dataframe(gdf)
lisa = esda.moran.Moran_Local(gdf['ALAND'], w_queen)

# ValueError: cannot assign slice from input of different size

I have shown here with EXAMPLE Dataset but actually I'm getting error with my original dataset also.
I did some research and couldn't find anything related to the problem. Moran_Local is working fine with other weights like KNN for my dataset but with Queen, it fails.


texts:
Is it possible to find out whether the user zooms in or out when using a Leaflet zoomstart or zoom event?

map.on('zoomstart ', function() {

    console.log('I am zooming, but am I zooming in or out?');

});


texts: Why do simple joins on small datasets slow down ArcMap so badly?

texts: ArcGIS 9.3.1, ArcInfo. I have some small file geodatabase and SDE featureclasses (few thousand records each, 50 columns) that I am doing attribute joins on. After the joins, drawing of the map is slowed down from sub-second to several seconds each time, sometimes more. Add in some definition queries that reference fields from the join and it's time to go get a cup of coffee while redrawing. Any ideas on why this happens? Anyone else ever experienced this?


texts:
I have a set of spatial objects in a MSSQL server 2008 R2 trying to dissolve adjoining objects where one of the attributes are equal. The table below is somewhat representative to my original table:

id     year     GEOM
01     1999     0x79F...
02     2013     0x79F...
03     1988     0x79F...
04     2013     0x79F...
05     1999     0x79F...
06     2013     0x79F...

Many of the objects above are neighbours so I would like to aggregate the spatial objects so that I in the end of the day end up with a smaller number of table entries/objects where possible. I am not looking for myltipolygons.

I have looked around and a propper dissolve functionality like the one available in for example ArcGIS desktop. MSSQL 2008 R2 does not seem to sport this functionality directly.

Is there a way to send a query to the database which does this job for me?


texts:
When I zoom to maximum; more than 18 the basemap is gone. How can we keep the basemap even when we zoom to maximum?
This is my plunker example


texts: Apply transparency to a raster layer by using actions

texts:
i have a specific action to load an raster layer by clicking on it, however i need to know if is possible to open that specific file and automatically apply a style.

My idea is when someone click and open that raster file, automatically he applies a transparency of 50%. Beeing "raster_path" an attribute field with URL directory.

I have the present code in Action (python):

import os pathin='[% "raster_path" %]'
layerID = str(os.path.splitext(os.path.split(pathin)[1])[0]) qgis.utils.iface.addRasterLayer(pathin,layerID)

i tried, with diferent results, saving a .qml file and apply to it, however for diferent raster show very diferent results with very bad resolution.

Is there a way to apply a transparency of 0.5 (50%) directly on the code above or do i need to set and field with .qml URL and load somehow on the same code?


texts: How to change a line Symbol layer type using PyQGIS?

texts:
enter image description here

I need to change a line style from Simple Line to Marker Line, and after that, use it as an arrow and change it's rotation and colors using PyQGIS.


texts: Fixing orientation in Net Topology Suite

texts:
I have a MultiPolygon that was created by reading from a GeoJSON string.

string geoJson = getGeoJson();
var rdr = new GeoJsonReader();
var multiPolygon = rdr.Read<MultiPolygon>(geoJson);

My problem - GeoJSON recommends but does not require that the outer ring be counter-clockwise and the inner rings be clockwise.

It looks to me as if GeoJsonReader, when handed a polygon with a clockwise outside ring constructs a MultiPolygon with a clockwise outer ring.

How can I detect, in Net Topology Suite, that a ring is clockwise or counter clockwise?


texts:
I am trying import shapefile with PgAdmin III during "QGIS geoalgorithms"

but I can't find mistake.

I don´t have "Shapefile to PostGIS Import Tool" and DB manager doesn´t work.


texts:
How can I set the Grass Region in Python without already having a map in the Mapset?
Basically I am looking for the equivalent to this window

Currently I am only able to extent an existing region, which leads in my case to very big regions.


texts: How to cache Geowebcache tiles from applications

texts:
I have enabled disk quota to cache tiles that produced by WMS. Geoserver admin panel Tile Caching->Tile Layers list of cached layers. I can preview here on a demo page.

http://localhost:8080/geoserver/gwc/demo/pgs:topas?gridSet=EPSG:4326&format=image/png

After preview I see that disk quate increased. This works fine.

I am using my layers on my OpenLayers application but used disk quata is not increasing.
application is using url: http://localhost:8080/geoserver/wms/.....
 Used tiles are not inserting to disk. Why does not this?

(I have enable Tile Caching->Caching Defaults->Enable direct integration with GeoServer WMS)


texts: How to address Geodatabase must be enabled as geodatabase admin user error?

texts:
I am trying to use "enable Enterprise Geodatabase" but I am getting the following error:

"Geodatabase must be enabled as geodatabase admin user.Failed to execute "

I am using postgresql 9.1.3 database.

I have created the database named geo using this command: createdb geo and sed user using: CREATE USER sde WITH SUPERUSER LOGIN PASSWORD 'sdepassword';  and then I connected to database as sde user.

How can I solve this problem?


texts: How to add layer to OpenLayers, but without it getting automatically added in LayerSwitcher Control?

texts:
I have a OpenLayers web map with LayerSwitcher control added. I need this to switch some of my layers. But I am also adding and removing some layers, which are temporary and I do not want them added to the LayerSwitcher Control.

So is there a way to add layer to OpenLayers, but without it getting automatically added in LayerSwitcher Control?


texts: How do I create centroids inside polygons in Python?

texts:
When creating centroids from polygons, ArcMap has an option to placing the centroids inside (within the bounds of) the polygons. I am attempting to create centroids in Python. I have tried:
df['geometry'].centroids

The problem is that some polygons are shaped crescent-like, or doughnuts as in the shown example so the resulting centroids are outside the polygons.
Is there a way of making the centroids fall inside their polygons in Python?

Using the function below suggested by @ThomasG77 does produce a good result.
df['geometry'].representative_point()

However, is there a way to influence the placement of this point so that it falls as centrally as possible within the bounds of the shape?


texts: Can I delete only tables, or only objects sharing a string of text, from a geodatabase?

texts:
I've got about 1000 tables in a geodatabse that I need to delete. They are named IntersectedBuffers_1, IntersectedBuffers_2, etc, and they're the only tables in the database.

Could someone please tell me how to delete these? I've been trying to do it from the python window. I can successfully delete individual tables with the following code:

import arcpy
arcpy.env.workspace = r"C:\Documents and Settings\Administrator\My Documents\ArcGIS\Default1.gdb"
arcpy.Delete_management("IntersectedBuffers_2")

But how can I delete all of them? I've tried things like

arcpy.Delete_management("IntersectedBuffers*.tbl")

... but I'm really just guessing! Is there a way to delete based on a shared text string, or based on teh fact that they're all tables?

Thanks in advance,

Jay


texts: Decimal Degree data erases itself between edit sessions in QGIS

texts:
I have a shapefile with point data.
I calculate the $x and $y of the points successfully

I save the changes and close the edit session.
Once I reopen the attribute table, the data after the decimal point is gone.

The field type is Decimal.

Why does this happen?


texts: PyQgs: unresolved reference for QgsProject() and QgsVectorLayer()

texts:
In my code I'm using the classes QgsProject() and QgsVectorLayer(). However PyCharm gives me the following error:

  Unresolved reference 'QgsProject' Inspection info: This inspection detects names that should resolve but don't.

The Project Interpreter is "Python3.7(firstplugin) C:\OSGeo4W64\bin\pyton-qgis.bat" which should be correct.

Any ideas how to solve this problem? I think I need to import something but I can't find what


texts: Loop every element in the same layer with a function using ArcPy?

texts:
I have the following code:

import arcpy
from arcpy.sa import *

arcpy.env.workspace="D:/..."
sp="StartPoint.shp"          #Dataset that defines the zones
zone_field="FID"             #Field that holds the values that define each zone.
value_raster="UDEM.tif"  #Raster that contains the values on which to calculate a statistic
x_shift = 1                  #increment of x coordinate of the point
TableName = "ZStT"             #Name of output table
arcpy.env.overwriteOutput = True
outTable = ZonalStatisticsAsTable(sp, zone_field, value_raster, TableName, "NODATA", "MEAN")
arcpy.AddField_management(TableName, "parameter", "FLOAT", 3, "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")
cursor=arcpy.da.SearchCursor(outTable, ["MEAN"])

for row in cursor:
    i = 0
    value = row[0]        #Create a cursor to query value in position of the start point
    while (value == 0):
        cursor2 = arcpy.da.UpdateCursor(sp, ["SHAPE@XY"]) #Move the point 1 meter to right
        for row in cursor2:
            cursor2.updateRow([[row[0][0]+(x_shift), row[0][1]]])
            i += 1      #"Moved step" plus one

        outTable = ZonalStatisticsAsTable(sp, zone_field, value_raster, TableName, "NODATA", "MEAN")
        cursor = arcpy.da.SearchCursor(outTable, ["MEAN"]) #query the raster value in the point's location
        for row in cursor:
            value = row[0]

    if value != 0:          #if(True)break the loop
       cursor3=arcpy.da.UpdateCursor(outTable, ["parameter"]) #parameter=i/value
       for row in cursor3:
           row[0] = (i/value)
           cursor3.updateRow(row)

The purpose of function is to iterate over coordinate change until the raster value in one single point's location meet the need (!=0),then put the parameter (moved step/raster value) in table output"ZStT".

Part of this question utilizes code in other question posted by me, but previous one accentuate how to change coordinate and this one focus on how to loop element. They're asking different stuff and I'm specifically told if there's new question post another one.

How to execute this function to element (point) of input layer separately since process to whole layers lead to change of all element simultaneously.


texts:
I just want to add combo box on OpenLayers map panel. I don't know how to add combo to the OpenLayers panel so currently I have set absolute location to combobox, like this http://jsfiddle.net/4zNH6/

<div id="map">
    <select id="combo">
    <option value="hyderabad">Hyderabad</option>
    <option value="pune">Pune</option>
    <option value="bangalore">Bangalore</option>
  </select>
</div>

#combo{
    position:absolute;
    left:10px;
    top:10px;
    z-index:10000;
    width:100px;
    height:20px;
    background-color:#FFFFFF;
}

Is this a right approach to place element on map??= or is it possible to add element in panel(OpenLayers.Control.Panel)?? and then add panel in map.


texts: Wrapping layer names or headings in ArcMap legend?

texts:
I've been struggling with this for a bit now but I cannot seem to find a way to wrap either layer names or headings in a legend. I can see how to wrap descriptions, labels, and make multiline legend titles but nothing for layer names and headings.

I have a long layer name/heading as pictured below and I cannot get it to run past the width of the symbols and labels. I can think of a couple of workarounds including moving this layer to the top and putting the layer name into the Legend title but would prefer not to have this item at the top of my legend.


texts: How to exit QGIS after running a python script?

texts:
I would like to run a python script via qgis --code myscript.py and then immediately exit.  I'm using iface.actionExit().trigger(), which kills QGIS when I run it from the python console, but not when put in the script that I pass to --code.

What's the right way to immediately exit?  I'm running QGIS 2.0.1

Update:  I've also tried sys.exit().  QGIS catches it and pops up a window with the following:

An error occured during execution of following code:
  execfile('myscript.py')

  Traceback (most recent call last):
    File "", line 1, in
    File "myscript.py", line 14, in
      sys.exit()
  SystemExit

Update: os.kill(os.getpid(), 9) works but it's a dirty hack and I'm looking for something better.


texts: Using Python to parse an XML containing GML tags

texts:
I am trying to write a python script which pulls GML tags out of an XML file and formats them into WKT for inserting into a PostGIS database. I have been successful in doing so for an XML containing a single part polygon using the follow code:

    rootElement = ET.parse("GMLExample_Polygon.xml").getroot()
    wkt = ""
    for subelement in rootElement.getiterator():
        for subsub in subelement:
            if subsub.tag == "{http://www.opengis.net/gml}X":
                x = subsub.text
            if subsub.tag == "{http://www.opengis.net/gml}Y":
                y = subsub.text
                point_for_pol = "%s %s, " % (x, y)
                wkt += point_for_pol
    wkt = wkt[:-2]

This code clearly won't work for multipolygons. I am unsure how to access the geometry for each polygon tag ("gml:Polygon srsName="BNG") separately and pull only the geometry nested under it. I am trying to use ElementTree, not sure if this is the best module to use? The XML is structured as follows:

<Order xsi:noNamespaceSchemaLocation="http://lalala.com/xml_polygon_order.xsd">
    <OrderRequest>
        <CustomerReference>999998</CustomerReference>
        <SiteAddress>
            <Premise>456</Premise>
            <Street>long street</Street>
            <Locality/>
            <Town>London</Town>
            <County>London</County>
            <PostCode>PN1 1PN</PostCode>
        </SiteAddress>
        <SiteGeography>
            <gml:Polygon srsName="BNG">
                <gml:outerBoundaryIs>
                    <gml:LinearRing>
                        <gml:coord>
                            <gml:X>452847.6009</gml:X>
                            <gml:Y>18596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                    </gml:LinearRing>
                </gml:outerBoundaryIs>
            </gml:Polygon>
            <gml:Polygon srsName="BNG">
                <gml:outerBoundaryIs>
                    <gml:LinearRing>
                        <gml:coord>
                            <gml:X>452847.6009</gml:X>
                            <gml:Y>18596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                    </gml:LinearRing>
                </gml:outerBoundaryIs>
            </gml:Polygon>
        </SiteGeography>
    </OrderRequest>
</Order>

Thanks for any help.


texts: Radio Propagation Modelling, from a GIS view point

texts:
I've had a request from another department involved in telecommunications to revamp elevation data into GTOPO30, DTED or SRTM formats as input into a program called Radio Mobile. They are working on adding infill sites to low signal areas. Radio Mobile looks to be a program with history and a quick peruse of the message archives and tutorial sites shows there is active user community.

The reliance on somewhat dated elevation formats, e.g. they don't seem to have heard of geotiff and do everything in geographic lat/long, got me wondering about the overlap between GIS and radio signal analysis. What communities and/or software systems are there out there approach Radio Propagation Modelling from a GIS point of view? What's your assessment of their GIS "savvyness"?

EDIT: This could also be phrased as, where does someone who knows a something about GIS but nothing about radio signal propagation go to learn more?


texts:
I am working on several layers that share similar symbology in ArcGIS for Desktop (Basic).  I'm trying to remove duplicate Symbology from the legend. Manually I can do it by opening the properties dialog box > selecting the symbology tab > selecting layers to remove (see jpg2), but this is time consuming.

This is ArcGIS for Desktop (Basic) 10.2.2 & python 2.7.5.

In jpg3+4 the attribute table for the 2 layers was added

Here the legend is built, but there is duplicate symbology as well and it is unecessary. Now, the maps contains 32 layers with a multiple duplicate symbology.

Finally, i want that the legend will composed of the two layers and will look like that:

I am looking for a solution both manually and via python script.


texts: Concatenate only when all fields to be concatenated have values

texts:
I have two fields that I want to concatenate.  Field1 contains null values and Field2 has values for every record.  I want to concatenate only when both fields have values.  I tried using a function with an if/else statement but that doesn't work as the resulting concatenation will return values from field2.

def concat_fields(field1, field2):
    if field1 == ""
        return

    else:
        return field1 + field2

I know this is wrong.  Any suggestions?

EDIT:


texts: St_ClosestPoint within St_MakeLine Query problem

texts:
I am currently trying to create the closest point (closest_ports_point_test) on the grid (lng_grid_test) from the original points (ports_point), and to create the line connecting the original and the closest points.

CREATE TABLE closest_ports_point_test
AS
SELECT *
FROM ST_MakeLine(
    ST_ClosestPoint(ports_point.geom, lng_grid_test.geom),
    ports_point.geom
    );

And I also tried :

CREATE TABLE closest_ports_point_test
AS
SELECT *
FROM ST_MakeLine(
    ST_ClosestPoint(
        lng_grid_test.geom::geometry(LineString, 4326),
        ports_point.geom::geometry(Point, 4326)
        ),
    ports_point.geom::geometry(Point, 4326)
    );

But I get the following error warning for both results:

  ERREUR:  entrée manquante de la clause FROM pour la table « lng_grid_test »
  LINE 2: AS SELECT * FROM ST_MakeLine(ST_ClosestPoint(lng_grid_test.g...
  ********** Error **********
  ERREUR: entrée manquante de la clause FROM pour la table « lng_grid_test »
  SQL state: 42P01
  Character: 310

Do you have any idea about what could be wrong with that?


texts: Converting/ translating Excel Formula in SQL query

texts:
i have a layer in QGIS "Network" with attribute table.  I want to compare columns rows from each columns and create ine and fill it based on result from comparison. Untill now I did it with help in Excel:

I'm looking for to this in QGIS with help of SQL query:

SELECT t1.*,

     CASE WHEN Rock_G_ID =
                `(SELECT Rock_G_ID`
                   ` FROM network t2`
               ` WHERE t2.rownum = t1.rownum - 1) THEN '1'`
     ELSE '-' END AS Class`
     FROM network t1`

This works for one column and can't get it for the both same time like the formula in Excel. How can I do this?


texts: Can GRASS v.net.iso split between nodes?

texts: Can GRASS's v.net.iso module split networks between nodes on the network, i.e., if I am finding the area within a cost of 400, and there is a line segment with a cost of 800 leaving the center node, will v.net.iso split that line halfway down it? (For the purposes of this question, the center node cost is 0).


texts: Embedded QGIS map functionalities in User Interface

texts:
I am using QGIS to analyze a set of data in which there are locations of different interest points and I need to calculate the distance between these locations considering roads network.

I have already done this adding the points as layers and calculating the shortest path with QNEAT3 tool.

Next I want to create a user interface so that user can add points to map, compute the shortest paths and then call an algorithm I have developed in C++.

Which is the best way to embedded QGIS functionalities namely interaction with map in a user interface? Could QGIS be integrated in a user interface made in C++ using Microsoft Foundation Classes?


texts:
I have a vector layer with lots of features and I'd like to create a new layer in Python, containing only a subset of those features.

To be more specific, I downloaded OpenStreetMaps data of Austria (http://download.geofabrik.de/europe/austria-latest.osm.pbf) and I'd like to extract features with "leisure"=>"playground" and put them to a new layer. I would like to create about 20 layers with various key-value pairs extracted from the original OSM data.

The problem is, that my current implementation is too slow:

layer = iface.layerTreeView().selectedLayers()[0]
layer_new = layer.materialize(QgsFeatureRequest().setFilterExpression('("other_tags" LIKE \'%"leisure"=>"playground"%\') OR ("leisure" = \'playground\')'))

On Austria's multipolygons layer, it takes about 2 minutes. On Corsica's multipolygon layer, which has much less features, it takes just few seconds. But I need it to run also on larger countries.

When I converted the Austria's OSM data to a SpatiaLite database and ran this query in the spatialite-gui client:

SELECT COUNT(*) FROM multipolygons WHERE other_tags LIKE '%"leisure"=>"playground"%' OR leisure='playground';

It took about four seconds (the result was 4708).

Is there any way to speed up the processing in the QGIS Python console or a script?


texts: Dynamic data on overlay layer in openlayers

texts: I am using openlayers to visualize an industry using image as the layer. On the raster layer I have created an overlay layer.Is it possible to show updated values (dynamic runtime values - as integer) in the overlay layer that gets updated every second? If so please help me how to get this done in openlayers.


texts: How to deploy a QGIS Python Standalone Application with minimum installation size?

texts:
I am trying to create a Standalone QGis Application using python. I was able to create the nice looking shapeviewer application by following the steps here.

Now I would like to deploy this in another system. Considering this is a small application, how can i keep the installation to a smaller size in the target machine? I would like to explore the possibility of not having to install QGis to run such a simple application. Is this possible?


texts:
I have a map with some layers, and the layers have different features on them. I have a function that iterates through each layer and creates an extent from all the features' geometry extent. Then I fit the map to this extent. I want to have a different center for the map and still fit all the features in the view of the map.

Is there a function that can enlarge/extend a given extent to have a new center?

I'm using OpenLayers 4.6.4
Example of code, below, so in essence I want to give a center parameter and have a function that extends/enlarges the extent already calculated.

 /**
 * Fits all the features in a layer on the map.
 * @param {*Layer on map which should get all the features from} layer
 */
function fitLayersFeatureOnMap(layers, center) {
    var featuresOnLayers = [];
    _.each(layers, function(layer) {
        featuresOnLayers.push(mapLayers[layer].getSource().getFeatures());
    });

    var extentOfAllFeatures = ol.extent.createEmpty();

    _.each(_.flatten(featuresOnLayers), function (feature) {
        ol.extent.extend(extentOfAllFeatures, feature.getGeometry().getExtent());
    });

    map.getView().fit(extentOfAllFeatures);
}


texts: Imported csv file from Mapinfo does not show any field related to color

texts:
I want to know how the color scheme works for imported file in Mapinfo.

I have a number of sites and links in Pathloss (they are colored in red and green for sites and links (e.g. On air and Not on air sites). I have exported them for Mapinfo. Mapinfo is showing the same color for sites and links as it was in Pathloss. It works good until now.

Now, I want to work with the thematic map (I need to show the legend of this layer as well as other layers) but can't find how can I choose the sites/links colorwise ( I have no column that shows on air/not on air or green/red or 0/1).

Is there anyway I can play with the data that is imported from Pathloss in thematic map?

Is there anyway to differentiate the sites in thematic map based on the color scheme they are using?


texts: QGIS: Get data in raster DTM from vectorpoints

texts: I have a vector layer of points. Can I get an output from a rasterlayer which contains DTM-values (height) determined by my points placement?


texts:
I am using eight day composite MODIS LST data MOD11A2. How do I make a monthly composite from four 8-day composite data taking the maximum value out of four in the output raster?
Initially I went for Raster Composite tool in ArcGIS 10.2.2 but that caused all void values to be 0K in the output which is not possible.
As I want values of all pixel for further analysis, how can I get the monthly composite for it?


texts: Exporting from ArcGIS Pro results in clipped output

texts:
I'm trying to export a world map from ArcGIS Pro 2.2.4, using the > Share > Map Export option. The map is zoomed to the full extent of the world:

But the resulting PNG or PDF is clipped to a square extent which doesn't cover the whole visible extent, or the extent of the datasets:

This occurs whether or not I choose to export a world file. This is my first time really using ArcGIS Pro so I'm just following the steps I would have used in ArcMap. Are there any steps I'm missing?

If it makes any difference I'm running MacOS with Windows 10 inside a Parallels virtual machine.


texts: How to implement workaround for ESRI BUG-000119907, arcpy.ApplySymbologyFromLayer_management

texts:
I'm trying to write an ArcPy script, within an ArcPy toolbox and the arcpy.ApplySymbologyFromLayer_management command is failing. as per bug report: https://support.esri.com/en/bugs/nimbus/QlVHLTAwMDExOTkwNw==

The workaround described however does not make sense to me.

  To fix this issue: 1. After adding the arcpy.ApplySymbologyFromLayer_management(layer, lyrfile) Python command, add the arcpy.SetParameterAsText(2, layer) parameter. 2. Save the script. 3. ArcGIS Pro > Catalog Pane > Toolbox > Right click Script > Properties. 4. Properties Window > Parameter Tab > Add a new parameter of type Layer, and set it to derived output

The layer is already set as a parameter of the script:

        param5 = arcpy.Parameter(
            displayName="TOP_10_LYR",
            name ="TOP_10_LYR",
            datatype="GPFeatureLayer",
            parameterType="Derived",
            direction="Output")
        param5.value = "TOP_10_LYR" ```
-------------------------
``` arcpy.ApplySymbologyFromLayer_management(TOP_10_LYR, TOP10_LAYER)```
------------------------------------------------
 arcpy.ApplySymbologyFromLayer_management(TOP_10_LYR, TOP10_LAYER) <-- not working!!

I have tried adding the line arcpy.SetParameterAsText(2, TOP_10_LYR)  after the ApplySymbology command but it doesn't seem to do anything.

Interestingly the first time I use ApplySymbologyFromLayer it does work, but not layers I create after that. Note that the layers are all new feature classes.


texts:
How can I force all labels to be visible in GeoServer

I have tried with both

<VendorOption name="conflictResolution">false</VendorOption>

and

<VendorOption name="spaceAround">-1</VendorOption>

But it doesn't seem to help much. The pointsymbolizer are always visible though.
And sometimes it doesn't even feel like its "crouded" when the labels isn't visible.

The (I guess) interesting parts in my sld file looks like following:

<PointSymbolizer>
    <Graphic>
                     <Mark>
                       <WellKnownName>circle</WellKnownName>
                       <Fill>
                         <CssParameter name="fill">#FFFFFF</CssParameter>
                       </Fill>
                      <sld:Stroke>
                        <sld:CssParameter name="stroke">#9F0000</sld:CssParameter>
                        <sld:CssParameter name="stroke-width">2</sld:CssParameter>
                      </sld:Stroke>
                     </Mark>
                     <Size>9</Size>
                   </Graphic>
                 </PointSymbolizer>
                <TextSymbolizer>
                <sld:Geometry>
                  <ogc:Function name="centroid">
                    <ogc:PropertyName>geom</ogc:PropertyName>
                  </ogc:Function>
                </sld:Geometry>
                    <Label>
                      <ogc:PropertyName>label</ogc:PropertyName>
                    </Label>
                    <Font>
                      <CssParameter name="font-family">Arial</CssParameter>
                      <CssParameter name="font-weight">Bold</CssParameter>
                      <CssParameter name="font-size">14</CssParameter>
                    </Font>
                    <LabelPlacement>
                      <PointPlacement>
                        <AnchorPoint>
                          <AnchorPointX>0.5</AnchorPointX>
                          <AnchorPointY>0.5</AnchorPointY>
                        </AnchorPoint>
                        <Displacement>
                          <DisplacementX>0</DisplacementX>
                          <DisplacementY>-15</DisplacementY>
                        </Displacement>
                      </PointPlacement>
                    </LabelPlacement>
                    <Halo>
                      <Radius>
                        <ogc:Literal>2</ogc:Literal>
                      </Radius>
                      <Fill>
                        <CssParameter name="fill">#9F0000</CssParameter>
                      </Fill>
                    </Halo>
                    <Fill>
                      <CssParameter name="fill">#FFFFFF</CssParameter>
                    </Fill>
                    <VendorOption name="conflictResolution">false</VendorOption>

                </TextSymbolizer>


texts:
I have a road network table. When I check st_isValid(the_geom)='f', I am getting an error:

  NOTICE:  IllegalArgumentException: point array must contain 0 or >1
  elements

  ERROR:  Exception in LWGEOM2GEOS

When I open this table in QGIS, I can view and edit the data. However, its attribute table is blank, while in PostGIS, the table has data.

What is the problem?


texts:
I am producing Kernal density estimates for animals from tracking data.

Location is input at x y data, I produce a spatial polygons dataframe from these coordinates.

xy<-id[c("X", "Y")]
id<-id[c("Individual")]
idsp <- data.frame(id)
coordinates(idsp) <- xy
proj4string(idsp)

I set the CRS

proj4string(idsp) <- CRS("+proj=utm +datum=WGS84 +units=m +no_defs")

I then use adehabitatHR to produce KDE95s for my individual animals using the kernelUD and getverticeshr.

cromKDE<- kernelUD(idsp, h = "href")
CromKDE95 <- getverticeshr(cromKDE, percent = 95)

I convert these outputs into the appropriate CRS for the package leaflet using

CromKDE95<- spTransform(CromKDE95, CRS('+init=epsg:4326'))

I then try to use the package leaflet to visualise the KDE95 polygons on a basemap

leaflet(CromKDE95) %>% addTiles() %>%
  addPolygons().

But the polygons are plotting in the Ocean south of San-Pedro rather than in Ireland as they should be.

I realise this is likely a simple fix adjusting the CRS but I have tried all the different formats I can think of.

Any help on this would be greatly appreciated.


texts: Ubuntu 'Manage plugins' in 1.9 (master) doesn't show 1.8 plugins

texts:
I'm using Qgis master on Windows 7 and more recently xubuntu 13.04. On xubuntu, the 'fetch python plugins' is missing, which Nathan just pointed out was intentional rather than a bug. My problem is that on xubuntu, Qgis searches the repository with unchangeable query switch

?qgis=2.0

This means the majority of plugins aren't available. Maybe I'm being really dumb here, but these work on my windows machine as I can install them through the python installer. Aren't 1.8 plugins mostly/all compatible with 1.9?
Is there any workaround?
When v2.0 hits, will all plugin developers have to have upgraded their plugins for compatibility?

Thanks!

Dez


texts: Import functions in python Add-In logic script

texts:
I am running into an error. I have a set of functions which I want to import in the python Add-In script at the time of execution, but when I built Add-In with "import stats" line then tools are unresponsive. Any suggestions whether it is possible to import functions in the Add-In script?

import arcpy
import os
import pythonaddins
import stats

optfolder = "C:/temp"

class AoI(object):
"""Implementation for rectangle_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.cursor = 5
        self.shape = 'Rectangle'
        os.makedirs(optfolder)

    def onRectangle(self, rectangle_geometry):
   """Occurs when the rectangle is drawn and the mouse button is released.
    The rectangle is a extent object."""

        extent = rectangle_geometry
        arcpy.Clip_management(r'C:/temp/ras', "%f %f %f %f" %(extent.XMin, extent.YMin, extent.XMax, extent.YMax),
                              optfolder + '/ras1', "#", "#", "NONE")
        arcpy.RefreshActiveView()


texts: I don't get any output from Reprojection

texts:
I have recently posted a Reprojection issue. I found the solution, but now I can't get the output: I load the resulting reprojection by code, but first, I look into the folder where I have indicated the path. There isn't any data ( I have tried "C:.." first)
I wonder if somebody could test this in order to confirm that this code works.. I think that I have a problem of system administration.

The code is as follows:

import processing
from processing.core.Processing import Processing
from qgis.core import *
from qgis.gui import *

Processing.initialize()
Processing.updateAlgsList()

rasterLayer = QgsRasterLayer("G:\..example...\SatImage.tif","MTN50")
print rasterLayer.isValid()
crs = QgsCoordinateReferenceSystem(25830,
       QgsCoordinateReferenceSystem.EpsgCrsId)
rasterLayer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayers([rasterLayer])#Loading the first     raster

processing.runalg("gdalogr:warpreproject", rasterLayer, "EPSG: 25830", "EPSG:4326", '0', 0, 0, 0, 0, 2, 1, 1, False, 1, False, "extraParam", "G:\..newFolder..\reprojected.tif")
reproject = QgsRasterLayer("G:\..newFolder..\reprojected.tif","ReprojLayer")
QgsMapLayerRegistry.instance().addMapLayers([reproject])


texts:
I tried osmosis, osm2pgsql, osmfilter, and probably most of OSM related tools in order to extract polygon rivers and water bodies.

Among the tools, OGR is intuitive and returns good results so far

Based on the OSM TagInfo, I am supposed to use the tags below to have river/water polygons

water=*, waterway=*

Here are the results of water=* from multipolygon below

It maybe look legit, though I do not see wide rivers (the basemap is OSM).

And I extracted waterway=riverbank from linestrings and it only gives me a partial riverbanks from Openstreetmap.

So, just to check what I've missed using query conditions, I extracted all the linestrings and checked the missing river areas.

It does not have any data as you may see.

I was just curious how come OSM has polygon shape of riverbanks/river polygons? if there are no information in raw data?
Or, is it OGR's failure?

And, are there any ways that I can have water areas from openstreetmap raw data?


texts: How can I load VBA customizations stored in a "document-level custom template" into a Map Document created from a "Normal.mxt" template in ArcGIS 10?

texts:
I have all kinds of custom tools that are called from UI Controls stored in a document-level map document (my custom template). In earlier versions of ArcGIS these were stored in the global-level "Normal.mxt". Since this it is no longer possible to save UI controls within the Normal template in version 10 I use the above mentioned map document that has all of my VBA customization to create most of my maps within so that they may "inherit" the customization. This is a very crude solution since it created a problem with maintaining and upgrading my custom tools until i migrate to .NET.

Is there a way to load document-level customization into a map document created from a Normal.mxt template?


texts: Calculate the average of neighbor pixels for raster edge

texts:
I have a small raster (or matrix) and I would like to know what is the difference between each pixel and the average of its' 8 neighbors in the case where the kernel window is on the edges and no values exist. When the 8 neighbor pixels have "real" values I can calculate the average of the 8 for each pixel, by using scipy.ndimage.convolve with kernel window as np.array([[0.125,0.125,0.125],[0.125,0,0.125],[0.125,0.125,0.125]]).

The problem is on the edges.

For example, if this is my matrix:

array([[13, 21, 13,  8],
       [ 5, 10, 22, 14],
       [21, 33,  9,  0],
       [ 0,  0,  0,  0]])

I would like to calculate the average of the neighbors of the upper left cell (value 13). The desired calculation will be

  (21+10+5)/3 = 12

For the pixel in the first column and second row (value 5) the desired calculation will be -

  (13+21+10+5)/4 = 12.25

Note that I need to ignored the central pixel (value 13 or 5) and all the cells that are outside of the 3x3 kernel window and also that the number of neighboring pixels with "real" values is different from pixel to pixel.

I read about the mode that deals with edges here and in other places but could not find a solution.


texts: Reclassification of environmental values

texts:
I've started to use Qgis recently.
I have a logical question.
I'm trying to reclassify a raster. The values in the raster are depths, obtained by an interpolation.
The original depth values (of the vector) are:
1 = <30°
2 = 30°-60°
3 = >60°
after the interpolation with IDW I get a set of continuous values and I need to convert them into discrete ones.
If I round up/down the continuous values I will obtain an over-estimation of the central values (2)
Is there a better (more logical) method to reclassify this values?


texts: Postgresql POSTGIS not using index field while having GIST index on geometry field

texts:
I have around 2M rows on my PostgreSQL table with GIST index on geo_point field but while I try to use ST_Distance or ST_DWithin it is not using the index instead it is doing Parallel Seq
Indexed Field:
CREATE INDEX rental_properties_geo_point_idx ON rental_properties USING gist(geo_point);

Initial Query: with ST_Distance
EXPLAIN SELECT id
FROM rental_properties
WHERE ST_Distance(ST_Transform(rental_properties.geo_point, 26986),
ST_Transform(ST_SetSRID(ST_MakePoint(-115.74181,46.394012), 4326), 26986)) / 1000 < 5
AND rental_properties.published = true
AND rental_properties.usd_price IS NOT NULL
AND rental_properties.usd_price != 0
ORDER BY distance ASC,
rental_properties.review_score DESC, rental_properties.number_of_review DESC
LIMIT 12

Explain:
->  Parallel Seq Scan on rental_properties  (cost=0.00..23188911.68 rows=197196 width=1016)
      Filter: (published AND (usd_price IS NOT NULL) AND (usd_price <> '0'::numeric) AND ((st_distance(st_transform(geo_point, 26986), '01010000206A6900004204EC0D715A47C1D41445CDFFE54041'::geometry) / '1000'::double precision) < '5'::double precision))

Then I converted it to ST_DWithin as per the suggestion here  like below but still, it is not using the Spatial Index
Converted Query: with ST_DWithin
EXPLAIN SELECT id
FROM rental_properties
WHERE
ST_DWithin(ST_Transform(rental_properties.geo_point, 26986),
           ST_Transform(ST_SetSRID(ST_MakePoint(-80.08713,26.237396), 4326), 26986), 5000::double precision)
AND rental_properties.published = true
AND rental_properties.usd_price IS NOT NULL
AND rental_properties.usd_price != 0
ORDER BY ST_DWithin(ST_Transform(rental_properties.geo_point, 26986),
           ST_Transform(ST_SetSRID(ST_MakePoint(-80.08713,26.237396), 4326), 26986), 5000::double precision) ASC,
rental_properties.review_score DESC, rental_properties.number_of_review DESC
LIMIT 12

Explain:
->  Parallel Seq Scan on rental_properties  (cost=0.00..17763977.78 rows=59 width=1009)
      Filter: (published AND (usd_price IS NOT NULL) AND (usd_price <> '0'::numeric) AND st_dwithin(st_transform(geo_point, 26986), '01010000206A690000899A6DDC3A0025C1AE1388EB44622AC1'::geometry, '5000'::double precision))

Another one suggested here to convert it from geometry to geography. So my basic question is how to do a radius search to utilize the index instead of Seq scan.


texts: Is there any method to retrieve SQL queries into a GeoJSON file using JSP?

texts: I want to display the contents of my PostGIS database of various shape-files in leaflet pop-up dynamically. I know there are static methods available that retrieves attributes. But I need to by passing SQL Queries to PostGIS database and retrieve the result as a GeoJSON file and display it in my map using leaflet. Is there any method to retrieve SQL queries into a GeoJSON file using JSP ?


texts: Selecting polygons below (downhill of) lines using ArcGIS for Desktop?

texts:
I am hoping to identify sections of a trail system which are close enough to potentially add sediments to wetlands downhill.

To this end, I would like to select polygons (wetlands) that are within a given distance (say 50m) below/downhill of polylines (trails) based on a 10m DEM raster.

Is there a way to buffer a line downhill a given distance? Or alternatively, is there a way to select by location in 3D? Any thoughts on potential workflows?

I am hoping that once I find a solution, I would be able to incorporate this process into model using ModelBuilder.

I am using ArcMap 10.2, have access to Spatial Analyst, 3D Analyst and most other tools.


texts: Printing map layers in pdf format using leaflet

texts:
I'm working with leaflet api to print map layers one by one. the output should be a pdf file with several pages, I want that in each page there will be a map (layer over a basemap)  on the right and a writing on the left. Is there a plugin that allows it with the ability to customize the output?
this is an exemple :


texts: QGIS mac: Problem running several algorithms from toolbox

texts:
I want to use the SAGA algorithm "add coordinates to points", but I get an error message:

  Algorithmus Add coordinates to points startet...
  shapes_points "Add Coordinates to points" -INPUT "/var/folders/ls/7vyvyj2x0pg4yrw2gtqgpfs80000gn/T/processing0f42e6b113ae40f2bb6eae1dac046156/1480657898.143.shp" -OUTPUT "/Users/Hirschi/Documents/Uni/Module/GEO113/Uebungen/Lab8/Data/coord.shp"
  /Users/Hirschi/.qgis2//processing/saga_batch_job.sh: line 3: 9194 Segmentation fault: 11 saga_cmd shapes_points "Add Coordinates to points" -INPUT "/var/folders/ls/7vyvyj2x0pg4yrw2gtqgpfs80000gn/T/processing0f42e6b113ae40f2bb6eae1dac046156/1480657898.143.shp" -OUTPUT "/Users/Hirschi/Documents/Uni/Module/GEO113/Uebungen/Lab8/Data/coord.shp"
  Converting outputs
  Lade Ergebnis Layer

  The following layers were not correctly generated.
  Points with coordinates
  You can check the log messages to find more information about the execution of the algorithm

The log message window looked like this:

  2016-12-02T06:51:38   0   SAGA execution commands
              shapes_points "Add Coordinates to points"  -INPUT "/var/folders/ls/7vyvyj2x0pg4yrw2gtqgpfs80000gn/T/processing0f42e6b113ae40f2bb6eae1dac046156/1480657898.143.shp" -OUTPUT "/Users/Hirschi/Documents/Uni/Module/GEO113/Uebungen/Lab8/Data/coord.shp"
  2016-12-02T06:51:38 0   Konsolenausgabe der SAGA-Ausführung
              /Users/Hirschi/.qgis2//processing/saga_batch_job.sh: line 3:  9194 Segmentation fault: 11  saga_cmd shapes_points "Add Coordinates to points" -INPUT "/var/folders/ls/7vyvyj2x0pg4yrw2gtqgpfs80000gn/T/processing0f42e6b113ae40f2bb6eae1dac046156/1480657898.143.shp" -OUTPUT "/Users/Hirschi/Documents/Uni/Module/GEO113/Uebungen/Lab8/Data/coord.shp"
  2016-12-02T06:51:38 2   Error loading result layer:
              Traceback (most recent call last):
                File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/gui/Postprocessing.py", line 75, in handleAlgorithmResults
                  out.name))
                File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/dataobjects.py", line 199, in load
                  + '\nCheck the processing framework log to look for errors')
              RuntimeError: Could not load layer: /Users/Hirschi/Documents/Uni/Module/GEO113/Uebungen/Lab8/Data/coord.shp
              Check the processing framework log to look for errors


texts: Creating Garmin GPS compatible map using TFW+tiff files

texts: I have complete set of 1:50000 maps of Sri Lanka in .tiff format. They have .TFW files accompanied with them. Now I need to create maps to use with my etrex 20 device using those tiff files. How can I do it?


texts: Find angle between two lines features that share a point

texts:
I'm trying to write a script that will find angles between polyline features. the polyline features always contain only two points (start and end).

A few sources indicates that the declination of a line can be calculated as so:

rads = math.atan2 (-dY, dX)

where dY is the difference in Y values and dX is the differece in X values of the start and end points. This seems to be where I'm failing.

Here's my code with comments:

import arcpy
import math

def GetAngle (p1, p2):
    x1, y1 = p1
    x2, y2 = p2
    dX = x2 = x1
    dY = y2 - y1
    rads = math.atan2 (-dY, dX) #wrong for finding angle/declination?
    return math.degrees (rads)

def LineToXYs (line): #return first and last coordinates
    firstX, firstY = (line.firstPoint.X, line.firstPoint.Y)
    lastX, lastY = (line.lastPoint.X, line.lastPoint.Y)
    return [(firstX, firstY), (lastX, lastY)]

def AngleFromLines (lines):
    #lines is a python list of line geometries that share a vertex
    for line1 in lines:
        for line2 in lines:
            if line1 == line2:
                continue
            line1StPnt, line1EndPnt = LineToXYs (line1) #get start and end xys for first line
            line2StPnt, line2EndPnt  = LineToXYs (line2) #get start and end xys for second line
            angle1 = GetAngle (line1StPnt, line1EndPnt) #calc angle - Doesn't work
            angle2 = GetAngle (line2StPnt, line2EndPnt) #calc angle - Doesn't work
            print "first line start and end coordinates:", line1StPnt, line1EndPnt
            print "second line start and end coordinates:", line2StPnt, line2EndPnt
            print "angle 1:", angle1
            print "angle 2:", angle2
            print "length 1:", line1.length
            print "length 2:", line2.length
            angle = abs (angle1 - angle2)
            print "angle between lines:", angle

AngleFromLines ([lineGeom1, lineGeom2]) #lineGeom* = arcpy polyline geometry object

For example, these two polyline geometries should return something around 270 degrees (or 90, also fine):

Instead I run the script and get:

first line start and end coordinates: (1974143.8040161133, 13562128.28540039) (1974120.358215332, 13562109.576782227)
second line start and end coordinates: (1974141.000427246, 13562131.794433594) (1974143.8040161133, 13562128.28540039)
angle 1: -0.000542988605953
angle 2: -0.000101843033066
length 1: 29.9952990964
length 2: 4.49148355856
angle between lines: 0.000441145572887

Any help would be great.


texts:
I am new to ArcPy and I am trying to write a code to select specific rows from a shapefile. I wrote this code, which gives back:

  Parsing error SyntaxError: invalid syntax (line 3)

Any suggestions about how to fix it?

    import arcpy
    arcpy.env.workspace = 'C:/data'
    arcpy.env.overwriteOutput = True
    fc = 'hospital.shp'
    with arcpy.da.UpdateCursor(fc, "Layer") as rows
        for row in rows:
            if row = "02-R00-04-B$0$210"
                pass
            elif row = "02-R16-04-B$0$210"
                pass
            elif row = "02-R00-04-B$0$220"
                pass
            elif row = "02-R16-04-B$0$220"
                pass
            elif row = "02-R00-04-B$0$211"
                pass
            elif row = "02-R16-04-B$0$211"
                pass
            elif row = "02-R00-04-B$0$320"
                pass
            elif row = "02-R16-04-B$0$320"
                pass
            else:
                rows.deleteRow()
    del rows
    del row


texts: How to correctly input complex multiband functions in Raster Calculator?

texts:
I'm trying to compute a MTVI2 (Modified Triangular Vegetation Index - Improved) and EVI (Enhanced Vegetation Index) from Landsat 8 imagery.

My formulas are as follow (I'm using ModelBuilder):

MTVI2=

1.5*(1.2*("%Band 5%"-"%Band 4%")-2.5*("%Band 4%"-"%Band 3%")) / SquareRoot(Square(2*"%Band 5%"+1)-(6*"%Band 5%"-5*SquareRoot("%Band 4%")-0.5))

EVI=

("%Band 5%"-"%Band 4%") / ("%Band 5%"+6*"%Band 4%"-7.5*"%Band 2%"+1)

I've tried threating the bands as floats with the same results. Any suggestions?


texts:
I am wondering how to join spatial polygons using R code?

I'm working with census data where certain areas change over time and I wish to join the polygons and the corresponding data and simply report on the joined areas. I am maintaining a list of polygons that have changes census to census and that I plan to merge. I'd like to use this list of area names as a lookup list to apply to census data from different years.

I'm wondering what R function to use to merge selected polygons and respective data. I have googled it but simply become confused by results.


texts: How to convert x y coordinates to OSM using QGIS?

texts: I would like to convert them to OSM file format (.osm).


texts: How can I import georeferenced photos (jpg, kml, kmz) using QGIS or ArcGIS for Desktop?

texts:
I've got several georeferenced digital photos as jpg, kml, and kmz.
How can I import them into my GIS projects (QGIS or ArcGIS for Desktop) so I can open them as a popup?


texts:
Intro

Short version: My goal is to find geometries intersected by a target polygon.

   Long version: We are building an interface for the user to upload new geometries (political maps) and allow him to solve 4D overlaps in a semi-automatic way.
 Steps: Find all intersecting polygons; Based on user review, modify polygons by using ST_Difference function.

To find overlaps - the most straightforward choice is to use ST_Overlaps or ST_Intersects. But this method will give us an unlimited amount of false-positive results. Probably because of the tolerance, borders located next to each other would be treated as an intersection.

SELECT *, ST_Intersects(
    ( SELECT territory FROM diff_test WHERE id = 1 ),
    ( territory )
) as intersects
FROM diff_test
WHERE diff_test.id != 1

The other method I can think of is to select every polygon limited by bounding box. But then there is a need to filter out all the "touching" borders, and ST_Touching is not helpful if you don't have common points. It will also not work.

Back to ST_Intersects:

Let's say we have two overlapped polygons - id 1 (1430) and id 2 (1487).

We've calculated two new polygons.

id 3 (1430 - 1487) ST_Difference(id1, id2)
id 4 (1487 - 1430) ST_Difference(id2, id1) # this one can be ignored

What did I expect:

id 3 is inside id 1 (ST_Overlaps = false)
id 3 is touching id 2 (ST_Touches = true)
id 3 does not intersect with id 2 (ST_Intersects = false)

Nothing from this was as I've expected.

The query was executed targeting id 3

Problem and Questions

What is the best solution to find truly intersect polygons?
What would be the best practice of avoiding intersections from newly created polygon by ST_Difference and his counterpart?

Thoughts

I can think of multiple ST_Difference operations, but it will not guarantee anything and it is slow. Next query still returns true instead of false

SELECT ST_Overlaps(t1mask, origt2) FROM ( SELECT
    (SELECT territory FROM diff_test WHERE id = 2) as origt2
    ,(SELECT ST_Difference(mask, intersects) FROM (
        SELECT
        ( SELECT territory FROM diff_test WHERE id = 3) as mask
        ,(
            SELECT ST_Intersection(
            ( SELECT territory FROM diff_test WHERE id = 3 )
            ,( SELECT territory FROM diff_test WHERE id = 2 ))
        ) as intersects
    ) as t ) as t1mask
) as t

As another option - cut a little bit more, to prevent this from happening.
Or maybe this problem can be reproduced only on my data and can be solved by ST_SnapToGrid or something similar

Test sample

Link for this shapes


texts:
I'm trying to georeference a GE satellite image and I've got the below problem. What I am doing wrong?

I want it as same as a square polygon to start mapping.


texts: How to find a vertex correct as from the latlong on pgRouting 2.1.0 +?

texts:
Using the syntax (SQL) then the 1º result found in my table of vertices is the street segment highlighted in blue, but if I use the "limit" greater than one (eg = 3), the street segment In red is found. The correct result would be the vertex of "Rua Expedicionário Guilherme Brasileiro", where is my lat long( current location). Observe the figure.

Using Limit = 1:

SELECT ways.gid,  ways.name, ST_AsText(the_geom) as the_geom FROM ways ORDER BY the_geom <-> ST_GeometryFromText('POINT(-34.935373 -8.037455)',4326) LIMIT 1;

Result:

65824, Rua São Mateus, LINESTRING(-34.9352902 -8.0367307,-34.9355167 -8.0373642)

Using Limit > 1:

SELECT ways.gid,  ways.name, ST_AsText(the_geom) as the_geom FROM ways ORDER BY the_geom <-> ST_GeometryFromText('POINT(-34.935373 -8.037455)',4326) LIMIT 3;

Result:

Can anyone suggest me other alternatives to find this information in my vertices table in a more precise way?


texts: Using "Endpoint Must Be Covered By" topology rule?

texts:
ArcMap 10.5.1
Windows 7 SP 1

Creating a drainage network. All runs of pipe are straight, and therefore must be capped on either end with a structure. So, I created a new GDB, FDS, and a set of topology rules.

Here's the issue: When I go to set up the topology (using the wizard in ArcCatalog), I am unable to select my point datasets from the list of "feature classes that will participate in the topology":


texts: How do I create blank geotiff with same spatial properties as existing geotiff?

texts:
How do I use gdal to create a blank geotiff with the same spatial properties (size, resolution, extent, projection) as an existing geotiff?

(The reason that I want to do this is that I'm using gdal_rasterize to rasterize a polygon shapefile to overlay with existing geotiffs, and I think it will be simpler to burn the values into a matching layer rather than specifying all the properties in the rasterize command.)


texts:
What would be the approach to show and dynamically update the map coordinates in the url like stamen is doing it? Are there any built in functions in OpenLayers, what other techniqes could I use or do I have to solve this on the server side?

Background: I have a simple TMS built with TileStache and a simple map for showing it - this all works fine. I was just wondering how to add this feature to my map


texts: Sourcing Layer in TOC to New Feature Class using ArcPy?

texts:
I'm trying to source a layer in my TOC with a new feature class. I've read countless forums, ArcGIS Resource Center, etc. With the MXD and python window open this is what I've surmised so far.

import arcpy

mxd = arcpy.mapping.MapDocumnet("CURRENT")

Using the following link, which scenario matches what I want to do most accurately?

MXD Source = r'\\na\data\telecom\projects\GIS\Proposal_Kentucky_Nextgen\02GISData\18DDP_ENG\MXDS\ENCHROACHMENT MXD\0180.SA - C001-C010 - 1.0 - PE.mxd

MXD Name = 0180.SA - C001-C010 - 1.0 - PE.mxd

Layer Name in TOC = STRUCTURE_SITEID_180

Layer in TOC's path = \\na\data\telecom\projects\GIS\Proposal_Kentucky_Nextgen\02GISData\18DDP_ENG\Design_Data.gdb\STRUCTURE_SITEID_180

Feature Class I Want to Source to Name = STRUCTURE_SITEID_211

Feature Class I Want to Source to File Path (same as the layer I'm trying to resource) = \\na\data\telecom\projects\GIS\Proposal_Kentucky_Nextgen\02GISData\18DDP_ENG\Design_Data.gdb\STRUCTURE_SITEID_211

How do I finish the Python script to accomplish the task? If I can get one to work I can extrapolate. I'm using ArcMap 10.3.1.


texts: How to update pgrouting stored procedure locally?

texts:
I think there is a bug in pgrouting 2.3 and would like fix it by changing my local stored procedure rather than waiting for the official release. I edited my local nodeNetwork.sql file, updated the extension, but did not see the changes; I also tried to drop and recreate pgrouting extension to no avail. What is the correct way to update and reload a PostgreSql extension?

As for the bug itself, I the following section in nodeNetwork.sql seems wrong:

--HAD TO CHANGE THIS QUERY
-- p_ret= 'create temp table inter_loc on commit drop as (
--        select l1id, l2id, ' || vst_line_locate_point || '(line,point) as locus from (
--        select DISTINCT l1id, l2id, line, (ST_DumpPoints(geom)).geom as point from intergeom) as foo
--        where ' || vst_line_locate_point || '(line,point)<>0 and ' || vst_line_locate_point || '(line,point)<>1)';
    p_ret= 'create temp table inter_loc on commit drop as ( select * from (
        (select l1id, l2id, ' || vst_line_locate_point || '(line,source) as locus from intergeom)
         union
        (select l1id, l2id, ' || vst_line_locate_point || '(line,target) as locus from intergeom)) as foo
        where locus<>0 and locus<>1)';

I believe the commented out part is correct, which was in the earlier release. The old code uses intersection point (geom) to cut lines, whereas the new code uses end points (source and target) of a line to cut anther line at intersection, which makes no sense. As the result, I see pgr_nodeNetwork() often creates 4 nodes for an intersection but not the intersection point, so no routes can be created between the 2 lines.


texts: pyqgis: using gdalogr:rastercalculator not working

texts:
I have an Int16 - 16 bit GeoTiff and would like to do a simple calculation.
Using the GUI in QGIS (raster calculator) it works without problems. The expression in the raster calculator is: Raster_Image * 0.04 - 273.15

I want to do the calculation via pyqgis. For this I would like to use the process gdalogr: rastercalculator.

I tried the following:

import processing
# Get actual raster layer (GeoTiff)
layer = iface.activeLayer();
output = processing.runalg('gdalogr:rastercalculator',
                                                layer,         #INPUT_A <ParameterRaster>
                                                '1',         #BAND_A <ParameterString>
                                                None,         #INPUT_B <ParameterRaster>
                                                '',         #BAND_B <ParameterString>
                                                None,        #INPUT_C <ParameterRaster>
                                                '',         #BAND_C <ParameterString>
                                                None,        #INPUT_D <ParameterRaster>
                                                '',         #BAND_D <ParameterString>
                                                None,        #INPUT_E <ParameterRaster>
                                                '',         #BAND_E <ParameterString>
                                                None,        #INPUT_F <ParameterRaster>
                                                '',         #BAND_F <ParameterString>
                                                'A*0.04-273.15', #FORMULA <ParameterString>
                                                '',          #NO_DATA <ParameterString>
                                                '1',           #RTYPE <ParameterSelection>
                                                '',         #EXTRA <ParameterString>
                                                None)        #OUTPUT <OutputRaster>
rLayer= QgsRasterLayer(output['OUTPUT'],'result')
print rLayer.isValid()

The calculated raster is not valid. I just do not understand what I'm doing wrong.


texts: Default value not updating on "Add Part"

texts:
I have a default value setup on a column in QGIS to calculate the area of a polygon. I also checked "Apply default value on update".
This works if I add a polygon or change the vertexes.

If I add a polygon part (Add Part) using advanced editing tools, the area is not updated.

Bug?


texts: Visualizing True Color Remotely Sensed Images

texts:
I am following this tutorial (which is also on GitHub) and I can't seem to understand this line of code:

np.clip(eopatch.data['BANDS'][closest_date_id][..., [2, 1, 0]] * 3.5, 0, 1)

I am guessing this is how you extract/retrieve RGB bands from a data-object, but I don't get this [closest_date_id][..., [2, 1, 0]] * 3.5, 0, 1 is it a standard thing? and what each component stands for?
This is what a single eopatch looks like, it's basically a stack of numpy arrays:

EOPatch(
  data: {
    BANDS: numpy.ndarray(shape=(28, 500, 500, 12), dtype=float32)
    CLP: numpy.ndarray(shape=(28, 500, 500, 1), dtype=uint8)
    NDVI: numpy.ndarray(shape=(28, 500, 500, 1), dtype=float32)
    NDWI: numpy.ndarray(shape=(28, 500, 500, 1), dtype=float32)
    NORM: numpy.ndarray(shape=(28, 500, 500, 1), dtype=float32)
  }
  mask: {
    CLM: numpy.ndarray(shape=(28, 500, 500, 1), dtype=uint8)
    IS_DATA: numpy.ndarray(shape=(28, 500, 500, 1), dtype=uint8)
    VALID_DATA: numpy.ndarray(shape=(28, 500, 500, 1), dtype=bool)
  }
  scalar: {}
  label: {}
  vector: {}
  data_timeless: {
    ARGMAX_B4: numpy.ndarray(shape=(500, 500, 1), dtype=int64)
    ARGMAX_NDVI: numpy.ndarray(shape=(500, 500, 1), dtype=int64)
    ARGMAX_NDVI_SLOPE: numpy.ndarray(shape=(500, 500, 1), dtype=uint8)
    ARGMIN_B4: numpy.ndarray(shape=(500, 500, 1), dtype=int64)
    ARGMIN_NDVI: numpy.ndarray(shape=(500, 500, 1), dtype=int64)
    ARGMIN_NDVI_SLOPE: numpy.ndarray(shape=(500, 500, 1), dtype=uint8)
    STF: numpy.ndarray(shape=(500, 500, 20), dtype=float32)
  }
  mask_timeless: {
    LULC: numpy.ndarray(shape=(500, 500, 1), dtype=uint8)
    VALID_COUNT: numpy.ndarray(shape=(500, 500, 1), dtype=int64)
  }
  scalar_timeless: {}
  label_timeless: {}
  vector_timeless: {}
  meta_info: {
    maxcc: 0.1
    service_type: 'processing'
    size_x: 500
    size_y: 500
    time_difference: datetime.timedelta(0, 7200)
    time_interval: ('2019-01-01T00:00:00', '2019-12-31T23:59:59')
  }
  bbox: BBox(((480000.0, 4050000.0), (485000.0, 4055000.0)), crs=CRS('32632'))
  timestamp: [datetime.datetime(2019, 1, 19, 10, 20, 45), ..., datetime.datetime(2019, 12, 25, 10, 20, 40)], length=28
)
More details on Eopatches can be found [here][2]

And this is the code for visualizing a single Eopatch

## Take a closer look at a patch

fig = plt.figure(figsize=(20,20))
eopatch = EOPatch.load(path_out +'/'+ f'eopatch_4', lazy_loading=True)
plt.imshow(np.clip(eopatch.data['BANDS'][closest_date_id][..., [2, 1, 0]] * 3.5, 0, 1))
plt.xticks([])
plt.yticks([])

fig.suptitle('Visible spectrum at single timepoint', y=0.93, size=20)
del eopatch


texts: Making Scripted Tools Load Output File Automatically into QGIS Project?

texts:
When writing your own scripts, and using some of the pre-fab tools how do you make the output load automatically into your QGIS project?

Example: In the Help Tab of the QGIS tool "saveselectedfeatures" it explains its "Console Usage" as follows:

processing.runalg('qgis:saveselectedfeatures', input_layer, output_layer)

But when you run it as its own tool from the Toolbox the GUI has you fill in the same two parameters, but also has a third: a checkbox for "Open output file after running algorithm".  I checked it and it showed up as promised on my map project.

But the Tool's own HELP tab says nothing about how to script this functionality.

So how do I make the output open in the project automatically?


texts: Geocoding Crime Dataset with Block Number only

texts:
I'm trying to geocode a list of crime statistics that include block-only addresses i.e.: 1900  BLOCK GRAND AVENUE
700  BLOCK REED AVENUE
3000  BLOCK MISSION BOULEVARD
and a couple thousand more for one zip code for one quarter period...

I have the city, state, zip as well.

I'm new to CARTODB and am working to add a layer of crimes to visualize along with other layers such as CENSUS tracts and another dataset of certain businesses. I have the zip, census tract and business layers added but need to add this layer of crime data.

I realize this has somewhat limited value because in a heat-map situation it would show the whole block for each individual crime. But, say I only pick out a few crimes like assaults. I would then have a low enough number where a heat map might make sense.

In any case the bigger question is I have to geo-code the block number for the address and am wondering if I can use street number ranges (i.e. 700-800) or a mid-block point (i.e. 750). The problem with both of these is the bias... i.e. the crime might have happened at the corner of 700.

I'm using the Texas A&M Geocode services (amazing service for free) https://geoservices.tamu.edu as well.


texts: Creating new rows from columns using FME

texts:
I am trying to convert columns to new features by using FME 2020.
I have tried Attribute splitter than feature joiner but until now without any great success.
The image below explains what I need.
Any tip ?


texts:
I'm trying to find the concave hull of a lat, lon pointset but when I project my points on my map it doesn't look like expected. I've tried a lineare projection and a spherical projection. My data is about 10000-70000 lat lon pairs. How is this possible?

$mapWidth    = 2000;
$mapHeight   = 2000;
$mapLonDelta = $mapLonRight - $mapLonLeft;
$mapLatDelta = $mapLatTop - $mapLatBottom;

$worldMapWidth=(($mapWidth/$mapLonDelta)*360)/(2*M_PI);
$f=min(max(sin($mapLatBottom*(M_PI/180)),-0.9999),0.9999);
$mapOffsetY=$worldMapWidth/2 * log((1+$f)/(1-$f));
$f=min(max(sin($mapLatTop*(M_PI/180)),-0.9999),0.9999);
$mapOffsetTopY=$worldMapWidth/2 * log((1+$f)/(1-$f));
$mapHeightNew=$mapOffsetTopY-$mapOffsetY;
$mapRatioH=$mapHeight/$mapHeightNew;
$newWidth=$mapWidth*($mapHeightNew/$mapHeight);
$mapRatioW=$mapWidth/$newWidth;

$tx = ($lon - $mapLonLeft) * ($newWidth/$mapLonDelta)*$mapRatioW;
$f = sin($lat*M_PI/180);
$ty = ($mapHeightNew-(($worldMapWidth/2 * log((1+$f)/(1-$f)))-$mapOffsetY));

This should look like Germany. What's wrong?:

Update: It seems to work. This is a map of Bonn(Germany):


texts: Difference between merge and create virtual raster

texts:
I have the DEM of a defined region that has been split into sub-DEM with a resolution of 1 km x 1 km. When I merge the sub-DEM to create a single DEM (Raster -> Miscellaneous -> Merge), the overall result is quite satisfactory, but for some reason, when I compare the original altitude to the merged altitude, there are differences in some regions. I also realized that the tool transforms regions containing nodata into 0.

The profiles in the following figure was obtained using the Profile plugin (red: original DEM, blue: merged DEM).

Is there an explanation for this?

I realized that when I create a virtual raster from DEM, and then save it in GeoTiff format, I get the original DEM.

My question is therefore: besides the ease of use, is there any advantage to using the Merge tool instead of the virtual raster?


texts: What does map scale refer to in map making / field survey metadata?

texts:
If I have, say, a printed map at the scale of 1:20000 and the legend says that it was mapped at a scale of 1:10000 does this simply mean, that the surveyors/cartographers made drafts & sketches or used aerial photographs at that scale or can there be some other meaning?

Also, if the metadata of vector data say that its scale/resolution is e.g. 1:5000, can it have any other meaning besides having been digitized from a source (aerial photograph, printed map) of that scale?

edit

See my comment below for specific context (general answers are still fine).