texts:
i was wondering

im getting different results when using geography() on different positions within a query.

So, whats the most accurate way of forcing geography() over a functions with geometry columns?

example:

SELECT

polys.id,

sum(st_length(st_intersection(geography(lines.the_geom),geography polys.the_geom))))/1000  total_lengths

from
schema.table1polys polys
schema.table2lines lines on ST_Intersects(geography(lines.the_geom),geography(polys.the_geom))

..or..

SELECT

polys.id,

sum(st_length(geography(st_intersection(lines.the_geom,polys.the_geom))))/1000  total_lengths

from
schema.table1polys polys
schema.table2lines lines on ST_Intersects(geography(lines.the_geom),geography(polys.the_geom))

whats the most accurate? i have sort of 2 / 3 Km of difference between both querys.


texts: Feature class SQL Query combining ANDs and ORs

texts:
I want to create a query that returns any combination of specific values of two different attributes. I have a feature class that has first field with values A through E and a second one that goes from 1 to 5. I want to select any features that are either A or B in the first field and either 1 or 2 in the second field. Would the following syntax be correct?

...WHERE ("field_1" = A OR "field_1" = B) AND ("field_2" = 1 OR "field_2" = 2)

or would it need to be making all the combinations?

...WHERE ("field_1" = A AND "field_2" = 1) OR ("field_1" = A AND "field_2" = 2) OR ("field_1" = B AND "field_2" = 1) OR ("field_1" = B AND "field_2" = 2)


texts: GDAL >= 2.2.x on Mac OS with Python bindings

texts: I'm needing to use the Python binding of GDAL >= 2.2.x in a Mac OS environment. Kyngchaos most recent binaries are provided for version 2.1.3. Out of compiling from source, is there any other option for install?


texts: Method selectedFeatures() works too slow in QGIS 3 Python console

texts:
I loaded spatial table(about 100 millions rows) from PostgreSQL(PostGIS) in QGIS3 and use method selectedFeatures() in QGIS 3 Python console. For such big tables it work very slow (it takes up to 10 minutes).

from qgis.core import *
from qgis.core import QgsVectorLayer, QgsDataSourceUri

uri1 = QgsDataSourceUri()
uri1.setUseEstimatedMetadata(True)
uri1.setConnection(str(host), str(port), str(database), str(user), str(password))
uri1.setDataSource (schema , tablename, geometrycol)
layer=QgsVectorLayer (uri1 .uri(False), tablename, "postgres")
QgsProject.instance().addMapLayer(layer)

iface = qgis.utils.iface
activeLayer = iface.activeLayer()
selected_features = activeLayer.selectedFeatures()

In QGIS 2 these methods work much more faster(several seconds).

How can I solve this problem? I want to get copy of selected data for several seconds.

Is this a bug or I am doing something wrong?

POSTGIS="2.5.1 r17027"

PostgreSQL 11.1 (Ubuntu 11.1-1.pgdg18.04+1)

QGIS version:3.8.3-Zanzibar

OS: Windows Server 2019

Python version: Python 3.7.0*


texts: What's the correct syntax for Execute SQL in Processing QGIS 2.6?

texts:
I made a small model with Processing in QGIS 2.6.

The final layer should be the result of a SQL query. Now, in the Execute SQL dialog window I can add the query:

select name, max(ele) as ele, geom
from layer_id
group by id

but I haven't understood which name have I to specify in the query (layer_id in the example).


texts: How to remove border line from one shapefile in interaction with another one

texts:
I'm new in QGIS, I have two shapefiles Northern Ireland (NI) and Ireland and there is overlapping with theirs borders.

Is there same way that could use border of NI for example and give me output merged those two files with border of one NI on inside part? Or can anyone help me how to solve this because as you can see I have also those "white parts" in between those two shapefiles.

I took advice from @Joseph, so I use shapefile (Ireland and NI together from ArcGIS and also the layer of NI constituency (it isn't county) NI_shapefile.zip. I did first

  Processing Toolbox > QGIS geoalgorithms > Vector general tools > Split vector layer

I got Group A (Ireland) and Group B (NI) from Ireland.shp
Then I use

  Processing Toolbox > QGIS geoalgorithms > Vector overlay tools > Intersection

using Group B (NI) and NI_shapefile (I change projection to be EPSG 4326) I get intersection layer. Now I wanted to check if everything is good and at first seemed like it is but then I zoomed and there is still white parts between borders  I tried to get layer with difference and then to do union but then QGIS crashes, every time.

I merged difference from NI_shapefile and Group B(NI that I got from splitting the whole Ireland) and added(merged) to Ireland part and that looks like this,

so I don't have white parts but I have like double lines in the boundary. Can I erase them with some function and use the outter boundary?


texts: am using Arcmap 10.0 to create a map. Need to display a chart for my data when users click on the marker. Must I use plugin like google chart? Can I do it directly in Arcmap?


texts: Display GeoTIFF layer from GeoServer in OpenLayers?

texts:
I have a number of GeoTIFF files to be published in GeoServer and consume in OpenLayers.

I have published vector data and I was able to consume them in OpenLayers. However, I am not able to show a sample GeoTIFF layer on OpenLayers (I created by placing the tiff file in GeoServer directory and published from there).

I used TileWMS source with Tile layer and ImageWMS source with Image layer, either way nothing gets rendered in OpenLayers. How can I display the GeoTIFF layer in OpenLayers?

code :

var url = 'http://localhost:8080/geoserver/localhost/wms'

var multipleSrc = new ol.source.ImageWMS({
    projection: 'EPSG:32643',
    url: url,
    params: {
        //'FORMAT': 'image/png',
        'LAYERS': 'localhost:Tiff',//getCommaSeperatedLayersFromHT(),
        //'TRANSPARENT': 'TRUE',
        //'CQL_FILTER': 'INCLUDE'//getCommaSeperatedLayersFromCQL() //"INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;tehsil='Jatara'" // Working
    },
    serverType: 'geoserver',
    //transition: 0,
    //transparent: true,
    //opacity: 0.5,
});

 var layers = [
    EsriLayer,
    //OSMLayer,
    multipleLayer,
    vector
  ];

  var view = new ol.View({
      center: washingtonWebMercator,
      zoom: 5,
      maxZoom: 19.9
  })
  var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: view
  });


texts: Lessen the overlapping in Polygons created on a road in PostgreSQL

texts:
I am working in an application of PostgreSQL. I have a list of points and want to create square of 500 sq. feet around these points. These points are in Geometry (EPSG:4326) format.Every point is center of square and final size of square should not be changed.

I have created a script in PostgreSQL to do this. But the problem is that on few places squares are more overlapped. I want to lessen the overlapping without changing there size and neither point should be outside of any square.Following data can be used for testing:

-84.26929 39.18387,-84.26929 39.18389,-84.26921 39.18399,-84.26613 39.18531,-84.26585 39.18542,-84.24252 39.19028,-84.24153 39.19059,-84.2408 39.19082,-84.23963 39.19123,-84.23803 39.19266,-84.23773 39.19336

Any idea?


texts:
When reading SRTM3 data from a .hgt file I get the wrong elevation of a specific latitude longitude and I can't see why.
Maybe someone can see the problem ?

My coordinate is near Lissabon in the sea and I expect an elevation around 0 or no elevation data, which means 0 as well. But I get an elevation of 32 meters. Something is slightly shifted, but I cant see what.

My coordinate :

My code in Java with values at the end of the lines for analysing:

Can someone see, what I am calculating wrong ?

Java code as text:

public double readElevation(LatLon coor) {
    System.out.println(coor);
    String tag = getHgtFileName(coor);

    ShortBuffer sb = cache.get(tag);

    if (sb == null) {
        return NO_ELEVATION;
    }
    // see http://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file
    double fLat = frac(coor.getLatitude()) * SECONDS_PER_MINUTE;
    double fLon = frac(coor.getLongitude()) * SECONDS_PER_MINUTE;

    // compute offset within HGT file
    int row = (int) Math.round(fLat * SECONDS_PER_MINUTE / HGT_RES);
    if(fLat < 0){
        row = row * -1;
    }else{
        row = HGT_ROW_LENGTH - row;
    }
    int col = (int) Math.round(fLon * SECONDS_PER_MINUTE / HGT_RES);

    int cell = (HGT_ROW_LENGTH * (row - 1)) + col;

    // valid position in buffer?
    if (cell < sb.limit()) {
        short ele = sb.get(cell);
        System.out.println("==> Read SRTM elevation data from row/col/cell " + row + "," + col + ", " + cell + " = " + ele);
        // check for data voids
        if (ele == HGT_VOID) {
            return NO_ELEVATION;
        } else {
            return ele;
        }
    } else {
        return NO_ELEVATION;
    }
}


texts:
Hey guys, I am looking for a real-time GIS vector engine that is to be used on a device using Embedded Linux/Qt (commercial or opensource).
Any suggestions? Thanks.

The usecase is a non-certified navigation system for airplanes, the data comes from a spatialite database. The data itself is not dynamic. Currently we are using rastercharts only, but we'd like to add the option of aviation vector charts, drawn on top of a raster layer.


texts: Splitting shapefile into separate files for each feature using QGIS or SAGA

texts:
I have a polyline shapefile with about 20 kinds of roads.
I would like to know how can I create automatically a vector file for each one.
I am aware of how to do it using ArcGIS Desktop - see Exporting feature class into multiple feature classes based on field values using ArcGIS Desktop?
How do I do it using QGIS or SAGA?


texts: Why are only certain features available for graduated symbology?

texts: I'm trying to use graduated symbology to show a color ramp for different ranges in the "Percent Ch" variable in my shapefile.  When I go to "graduated" in the layer styling panel it shows me that only some of the variables in my attribute table are available from the drop-down menu.  "Percent Ch" is listed in the attribute table, but I can't use it for graduated symbology.  The only thing I can see that's different with this variable is that it consists of whole numbers with decimals, while the other variables are only whole numbers.  Do you know why I can't use graduated symbology here? Do I have to convert the data into a different format?


texts:
I have thousands of GeoTIFF source files that cover whole country in 6 different layers.
I have created VRT file for each layer.
I loaded everything in QGIS, assigning correct colors to bands in each VRT.
I exported from QGIS a VRT of it all.
When running the combined VRT through gdal2tiles, the resulting tiles are all either completely blue, or part lightblue. I suspect this corresponds to two different original sources of layers.

Is there a way to preserve the color palettes as defined in QGIS? The challenge is, that each layer is composed of hundreds of monochrome files, and it seems the color assignment wasn't saved?


texts: Need to convert points to polygons of a specific size, using the point as the center of the polygon

texts:
I have a large collection of aerial photos with no georeferencing other than the Lat/Long location of the camera. Georeferencing them invididually is not an option (too much time, too many photos, some are oblique, some vertical). I need to automatically generate some kind of approximate polygonal footprint of each photo. I was thinking there must be some way to auto-generate a rectangle from a point layer using the point as the center and using numbers in the attribute fields as the dimensions. Or, barring that, just generating a circle at each point using an attribute as the diameter.

I'm using Arcview 9.3.1.


texts: Measuring polyline lengths running through a raster - is this possible?

texts:
I have a reclassified raster layer with 4 values and a vector layer of polylines. I would like to discover the length of the lines running through each raster value, or the percentage of the lines within each raster value.

Can anyone tell me if this is possible please?

I was informed Zonal Statistics may be able to do this, but to no avail.

Thanks.


texts: How to auto-fill a new point vector in a vector layer?

texts: When I work in-situ I manually add points in my shape-file. I would like that when I add a new point some attributes would already filled with values calculated programmatically. E.g in my first column a progressive number (first point automatically add 1, second 2 etc. etc.) like an auto $rownum second column id of the instrument measure usually is something like $rownum+nwhere nis a number. I'm able to to that only in post processing with the command I've mentioned but for me would be very helpful have a quick look to labeled point when I'm taking them.So is there a way to pre-fill determinate values of a new point vector?


texts:
I have a problem with displaying an image with QGIS2Web. First, I followed the procedure as given in the link below to pop some images from my photo directory:

Adding hyperlinks to pop-up data in qgis2web

Then I followed the step and typed like this:

  CONCAT('<img src = "E:/haris pictures/TL0101.JPG'  width="300" height="225"/>')

I put that code in my attribute table (FOTO HTML).

After that, I used QGIS2Web plugin to produce webpage for my map and it doesn’t end up with what I expected, like this:

Does anyone have suggestions to solve this problem?


texts: Increment VALUE column in raster in ArcGIS using ArcPy

texts:
I would like to increment the VALUE field in a raster by a constant amount (say 100).

I can do this by using reclassification, but wondering if there is a simpler more direct method (using ArcPy)?


texts: Python to change file names in a directory

texts:
I am converting a directory of kml files to 1 single fGDB, but they all have the date as the file name.  This creates a problem when using the original name because of the number at the beginning of the file name, it also has hyphens in the name of the file i.e. 2011-09-16 123.kml

I need to use python to change the - to _ and add KML to the beginning of every file in the directory.


texts: I want to add data in my map from ArcGIS Online, but ArcGIS does not allow me to complete that operation. I tried on the Add In Manager, by using the option Load All Ads-In without restriction, but it still does not allow me to add data from ArcGIS Online, nor from Add Basemap. I also tried to connect ArcGIS  to a ArcGIS Server, but I can´t find any server address in order to complete that operation. Any suggestions on how I can address this problem, in order to be able to Add Basemaps and Maps from ArcGIS Online?


texts: Finding angle between intersecting features in two feature classes using ArcGIS Desktop and Python?

texts:
I have two intersecting line featureclasses. I want to find the angle at each point of intersection using ArcGIS 10 and Python.

Can anyone help?


texts: I tried creating a space time cube using "create space time cube from defined locations" using a 1 day time step interval just to create a 3D visual. it runs successfully but the output has a creates cubes with values of zero for the dates of "missing data" the data set is composed of data across three years but only of January, February and parts of March. Is there a way to set this up so that I don't have to query out the zero data?


texts: Converting x,y coordinates (EPSG:102002) to latitude,longitude (EPSG:4326)

texts:
I have ranges of x and y coordinates in EPSG:102002 GRS 80 that I need to convert into EPSG:4326 WGS 84. I have tried different suggestions here and here but nothing has worked so far.
I find myself manually converting points of x and y coordinates to lat and lon coordinates online, one point at a time (N = 8200).
Is there a way to automatize this process where it also takes the projection into consideration?
I am pretty flexible with the program using, as long as it is not ArcGIS and Java.
Here is a sample
x             y            lat       lon
-3145588.03  5018975.86   63.30     177.62
-3115010.95  5010633.39   63.55     177.94
-3084474.59  5002315.52   63.80     178.28


texts: How to use if statement to check if a tag is assigned to particular value

texts:
I would like to filter out areas in Germany with post code = 14480. I coded the lines posted below
but when I run it returns nothing.
How can I fix this?
Code:
[timeout:3600];
area["ISO3166-1"="DE"][admin_level=2]->.germany;
foreach .germany -> .reg (
  node.germany(if: t["addr:postcode"] == reg.u(t["addr:postcode"]))
  ->.targetArea;
  node.targetArea(if: node.targetArea["addr:postcode"] == "14480");
  out geom;

);
  out;


texts: Building Spatialite network with long integer node?

texts:
After I build a spatialite network, I found that some node can't be routed. And after some survey, I found when the NodeFrom and NodeTo is bigger than 2^31 - 1, it will return NULL.

Is there any solution I can build the network with large integer?


texts: QGIS 3 how to customize line dash symbology

texts:
I want to keep the same line width but specify a much longer dash.

The symbology controls in QGIS 3 have changed significantly from 2.X, so the many answers re: line symbology control for QGIS 2.X are now obsolete.

When I attempt to edit the settings using the button-dropdown to the right of the Stroke style field, I invoke an Expression String Builder that might enable me to customize stroke styles to my heart's delight; but the syntax is not explained.

It looks like it expecting an SQL expression, because it looks like a similar dialog box for managing Attributes. But I cannot find any documentation on how it works. The "Expected Format: string [no|solid|dash|dot|dash dot|dash dot dot]" is cryptic.


texts: Can't install pyspatialite with python3 and sqlite3

texts:
I can't install pyspatialite. When I try it gives me this error:

pip install pyspatialite
Downloading/unpacking pyspatialite
Could not find a version that satisfies the requirement pyspatialite (from versions:
2.6.2-spatialite.2.3.1, 2.6.2-spatialite.2.4.0-4, 3.0.1-alpha-0)
Cleaning up...
No distributions matching the version for pyspatialite
Storing debug log for failure in /home/jenia/.pip/pip.log

The python interpreter that's installed in the virtualenv is 3.4.1.

When I run explicetly:

pip install pyspatialite==2.6.2-spatialite.2.4.0-4

I get this:

Running setup.py (path:/home/jenia/pip/app/build/pyspatialite/setup.py) egg_info for package pyspatialite
Traceback (most recent call last):
  File "<string>", line 17, in <module>
  File "/home/jenia/pip/app/build/pyspatialite/setup.py", line 83
    print "Is sphinx installed? If not, try 'sudo easy_install sphinx'."
                                                                       ^
SyntaxError: invalid syntax
Complete output from command python setup.py egg_info:
Traceback (most recent call last):

File "<string>", line 17, in <module>

File "/home/jenia/pip/app/build/pyspatialite/setup.py", line 83

print "Is sphinx installed? If not, try 'sudo easy_install sphinx'."

                                                                   ^

 SyntaxError: invalid syntax

 ----------------------------------------
 Cleaning up...
 Command python setup.py egg_info failed with error code 1 in /home/jenia/pip/app/build/pyspatialite
 Storing debug log for failure in /home/jenia/.pip/pip.log

Does anyone know how to fix this error?


texts:
I would like to create random points across my study area (which is a DEM raster) using the Feature Class> Create Random Points. However, I need to select/copy/move/delete these points and re-arrange them in a line, with each point 1km away from each other. I was wondering whether there is a specific tool that I can use in ArcGIS 10.1 to achieve this -I mean a tool that can tell me how far each point is from another point or is there way to measure the distance between two points to ensure I am placing these points exactly 1km away.

Any info would be appreciated!


texts:
I am trying to validate the geometry of all polygons in my table.

When I perform the following query, I get an error:

update parcel2016
set geom = St_MakeValid(geom)

ERROR:  Column has M dimension but geometry does not
********** Error **********

ERROR: Column has M dimension but geometry does not
SQL state: 22023

How might this be corrected? I've also tried wrapping in:

set geom = st_multi(st_collectionextract(st_makevalid(geom)))

but this was equally unsuccessful


texts: Assign “parcel boundaries” to “census areas” from “MULTIPOLYGON” in GeoPackage file “.gpkg” using Python

texts:
In one layer I have “polygons layer” of all “parcel boundaries”, in another layer I have the “polygon” of “census areas”. And I would need to assign every “parcel boundaries” to the “census areas” belongs to. I have a GeoPackage file “.gpkg”. Geometry is in MULTIPOLYGON format. E.g. MULTIPOLYGON (((595568.124 4417209.48, 595563.754 4417209.48, and so on…..)))
Do you know how can assign parcel boundaries to census areas from MULTIPOLYGON format using Python or find more information?
I am reading “.gpkg” file using “geopandas” in Python.


texts: Select for a specific number of characters

texts:
I am using ArcGIS 10.3.1 and have an attribute that has either 4 or 5 characters; i.e., CO-1, BY-3, CAR-1, MRT-8.

All letters are uppercase, this is in a personal GDB, and this is a string field.

I am trying to select for only the attributes that have 4 characters.


texts: I have a polygon grid file and want to clip my raster imagery to the shapes outlined by it. I am using arcmap and was wondering if there was a tool directly for this? I would like to also be able to name each resulting raster for the polygon feature it was clipped to.


texts:
This seems like it should be straight forward, but I'm struggling to get the correct transformation from my state plane to lat long.

My code is as follows

# projection from spatialreference.org/ref/epsg/2992/
# oregon is espg ref 2992
projection = r'+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048 +no_defs'

osp2992 = Proj(projection)

Now when I try and check against the input/output they show in the upper right corner

# Input Coordinates: -119.88481361812, 43.453242550659
# Output Coordinates: 1475673.108308, 621492.180223

>>>osp2992(-119.88481361812, 43.453242550659)
(449785.16341222316, 189430.81653184193)

It's way off.

I'm not sure where I went wrong.

Edit: removed the to_meter conversion and still getting wrong output

projection = r'+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 +datum=NAD83 +no_defs'

>>>osp2992(-119.88481361812, 43.453242550659)
(550261.2548834425, 157725.95224630416)


texts:
I want to add this extension to an sddraft using minidom and can't figure out how to do it.
Some help would be great.
To start with I'd like to find out how to add the first node SVCExtension xsi:type="typens:SVCExtension" right after Extensions xsi:type="typens:ArrayOfSVCExtension".

<SVCExtension xsi:type="typens:SVCExtension">
  <Enabled>true</Enabled>
  <Info xsi:type="typens:PropertySet">
    <PropertyArray xsi:type="typens:ArrayOfPropertySetProperty">
      <PropertySetProperty xsi:type="typens:PropertySetProperty">
        <Key>supportsSOAP</Key>
        <Value xsi:type="xs:string">false</Value>
      </PropertySetProperty>
      <PropertySetProperty xsi:type="typens:PropertySetProperty">
        <Key>supportsREST</Key>
        <Value xsi:type="xs:string">true</Value>
      </PropertySetProperty>
    </PropertyArray>
  </Info>
  <Props xsi:type="typens:PropertySet">
    <PropertyArray xsi:type="typens:ArrayOfPropertySetProperty"/>
  </Props>
  <TypeName>StreetDirectory</TypeName>
</SVCExtension>


texts:
I'm new to GIS but have been learning QGIS for the last few months. I have Windows10/64bit system and I just installed the latest QGIS version 2.18.10, but now I get an error when attempting to run GRASS plugins from the processing toolbox. Error is as follows:
"This algorithm cannot be run :-(
It seems that GRASS GIS 7 is not correctly installed and configured in your system. Please install it before running GRASS GIS 7 algorithms."

I have seen many post on this topic and have tried every suggestion to fix it to no avail.

I have completely uninstalled, including cleaning up the registry, re-installed 2.18.10 and still same issue.

I'm about to uninstall again and revert back to 2.18.08 which was working prior to this last update. But want to reach out and see if someone has a better solution prior to doing that.


texts: Downsampling Geotiff using summation - Gdal/Numpy

texts:
I am trying to downsample a 1 square km raster dataset to a much larger (.5 degree x .66 degree) dataset by summing all of the pixel values within this large grid cell.

gdal_warp does not contain a summation resampling method so I'm wondering if anyone has figured this out before.


texts: Debug QGIS3 plugin with Visual Studio Code

texts:
I would like debug QGIS plugin using Visual Studio Code. I have Win10 and QGIS 3.4.4. I followed the instructions written in these two sites:

How to setup QGIS 3 and Visual Studio Code plugin development (Windows)
Debug Visual Studio plugin homepage

When I start the debug in visual studio code I have the follow error message:

Exception in thread ptvsd.Server:

  Traceback (most recent call last): File "C:\PROGRA~1\QGIS
  3.4\apps\Python37\lib\site-packages\ptvsd\ipcjson.py", line 269, in process_one_message  msg = self.__message.pop(0) IndexError: pop from
  empty list

During handling of the above exception, another exception occurred:

  Traceback (most recent call last):
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\threading.py", line 917, in _bootstrap_inner
          self.run()
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\threading.py", line 865, in run
          self._target(*self._args, **self._kwargs)
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\wrapper.py", line 712, in process_messages
          raise exc
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\wrapper.py", line 697, in process_messages
          self.process_messages()
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\ipcjson.py", line 258, in process_messages
          self.process_one_message()
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\ipcjson.py", line 272, in process_one_message
          self._wait_for_message()
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\ipcjson.py", line 154, in _wait_for_message
          line = self._buffered_read_line_as_ascii()
        File "C:\PROGRA~1\QGIS 3.4\apps\Python37\lib\site-packages\ptvsd\ipcjson.py", line 113, in _buffered_read_line_as_ascii
          temp = self.__socket.recv(1024)
      socket.timeout: timed out

A similar error was reported in ptvsd github site, the bug was resolved, but I have already the error.
Did anyone have the same issue?


texts: Irregular temporal data and visualization

texts:
I want to add timed data to ArcGIS. The problem is that the data is very bulky and in Excel file there are two columns of date and time of start and end of the interval and these intervals are irregular, for example some under one hour, some are hours and some are days.
Is there a way to automatically introduce these two columns to the software and software identify start and end by itself?
Is it possible to show these data on time slider while the intervals are not equal?


texts: Migrate MySQL spatial DB/tables to PostgreSQL

texts:
We have a web system that uses MySQL spatial databases. We want to migrate from MySQL to PostgreSQL. I've been looking for methods to do it. I have found some tools that help you to migrate. However most of them do not support spatial data like https://github.com/philipsoutham/py-mysql2pgsql.
I also took a look over How to migrate spatial tables from MSSQL to PostGIS. I have just seen this post. I haven't tried it yet, I don't know if it could work for MySQL.


texts: Adding compression variable to CopyRaster_management for TIFF Raster Dataset?

texts:
I am using ArcGIS Desktop 10.1.

I need tiff orthos to have NO compression. I'm using a python script to parse through directories and add 4-band data to the images. It works great ... but the customer wants it uncompressed. The process creates them in LZW by default.

The copyraster info page states: out_rasterdataset --> "When storing a raster dataset to a JPEG file, JPEG 2000 file, TIFF file, or geodatabase, you can specify a compression type and compression quality."

How do I set the compression type to None?

I have tried arcpy.env.compression = "NONE" but this does NOT work. I even added the pyramid environment variable. No dice.

Here's what I'm using now and it works, but the compression is still LZW:

import os,arcpy,glob

rasterpath = r"source-directory"
outFolder = r"end-directory"

arcpy.env.workspace = rasterpath
arcpy.env.pyramid = "PYRAMIDS -1 BILNEAR NONE NO_SKIP"
arcpy.env.compression = "NONE"

for rasterFile in arcpy.ListRasters("*.tif"):
    oName, oExt = os.path.splitext(rasterFile)
    outRaster = os.path.join(outFolder, oName+ ".tif")
    arcpy.CopyRaster_management(rasterFile,outRaster,"DEFAULTS","255","255","","","8_BIT_UNSIGNED")

In order to get no compression, I have to open each file and export data. If the no compression variable won't work in the current script, I'm fine with creating a new script. Basically, I want this screen in script form:


texts: How to get link of wms layer in Geoserver

texts:
How to Get WMS links of all published layers from specific workspace and store in Geoserver using gsconfig  or any environment.
I want to get all WMS links of all layers in excel? I can get a published layer name from below code.

from geoserver.catalog import Catalog
cat=Catalog("http:// localhost:8080/geoserver/rest")
myResource = cat.get_resources(store='storeName', workspace='workspaceName')
layers = cat.get_layers(myResource)
for rss in myResorce:
  print rss.name


texts: Does a value in srtm dem file refer to the center or corner of a cell?

texts:
I'm working with srtm dem files. They are ascii text files containing a table/grid of elevation values. At the top of the file is some metadata to locate the grid in lat/lon:

xllcorner     -124.00041606132

yllcorner     48.999583599682

cellsize      0.00083333333333333

So initially I thought of the cells as points and assumed that the first value in the grid was the elevation for a point located at (xllcorner, yllcorner), and that the second value in the grid was the elevation at (xllcorner + cellsize, yllcorner). I then used this data to draw a 3d mesh with the mesh vertices corresponding to cells in the grid.

The resulting mesh looks correct but now I'm wondering if it's out by (cellsize / 2) in lat and lon...

When I started working with the file using Postgis Raster I discovered that it considers each cell as a polygon rather than a point. Rather than the first value representing the elevation for the point (xllcorner, yllcorner) it actually represents the the elevation for the polygon (xllcorner, yllcorner, xllcorner + cellsize, yllcorner + cellsize).

So, with that in mind, and seeing that I need to map each elevation to a point. Is it more accurate to say that the first value in the grid represents the elevation at (xllcorner, yllcorner) or at (xllcorner + cellsize / 2, yllcorner + cellsize / 2)?


texts:
I am trying to work with cordex climatic data (curvelinear and rotated data) in the stars package, and I am having a lot of issues, so decided to simply reproduce online example, and I fail reading in data with the command read_stars.

I am trying reproduce in this post https://github.com/r-spatial/stars/issues/52#issuecomment-430753055 by using the provided data here https://github.com/r-spatial/stars/issues/52#issuecomment-430548580

library(stars)
f = "1-RawData/example_ob_tran.nc" #path to the downloaded file

(r = read_stars(f))

When I do this I am getting

  Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver), : file not found

while the file the following command works

(r = read_ncdf(f))

suggesting that the file exist.

This seems very basic, but I am unable to sort this myself. Does anyone see what I am doing wrong?


texts: Conditional Con tool setting rgb value for true statement

texts:
I have a raster with some NoData values. I am trying to set these values to black. I have looked at the documentation.

This is what I have so far

Con(IsNull('MyImage'),0,'MyImage')

This produces a grey scale image. I want to maintain the color in the image and just set the NoData values to black.

I am using ArcMap 10.3


texts:
I have a gigantic OSM geodata file of a city, and I want to split it into smaller files according to each group of  tags.

I googled it up and i found that there is a jar file caller splitter.jar . i downloaded it and made a test with a small osm geodata, but what i got is a .pbf file which is not readable.

is there any recommended software to split the osm data or i have to do it programmatically?


texts: I'm building a map that uses a conditional statement to determine what data to load based on the zoom level. However, the conditional statement is only fired when the map first loads, as opposed to every time the zoom level changes. Is there a way to change the code to run on zoom load (i.e. to update every time the user zooms in or out) rather than the initial page load?


texts: How to download images from google maps using RgoogleMaps

texts:
I'm trying to download images from Google Maps using the package "RgoogleMaps" in R (getmap function). But I need to download a lot of images, so I've been thinking to create a fishnet in QGIS and import the coordinates to R, but I don't know how to import it to create a loop and download the images.

Does anyone know how to do this?


texts: QGIS Way to join based on different length keys fields?

texts:
In QGIS is there a way to join a .shp file to a table when the keys are different length (by ignoring part of one set of keys?)

For example .shp file key col like this
1500000US120730002001
1500000US120730002002
1500000US120730002003
(the 1500000US may be ignored because in this .shp file they are all the same)

Join File has a key col like this
120730002001
120730002002
120730002003


texts: Projection WGS 84 is slightly displaced referring to a hillshade in QGIS

texts:
I have a hillshade and a vector layer. The vector layer is slightly displaced referring the hillshade. About lat=100m and long=175m

The CRS of the project is EPSG:3397 - PD/83 / Gauss-Kruger zone 4.

The CRS from the vector layer: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs.

How can I place the vector layer on the right place in QGIS
I tried it with the +towgs84=... attribute. But this didn't work.


texts: Sentinel-1 image decible in Google Earth Engine?

texts:
I've wrote this code for Sentinel-1 decible image. although it doesn't show any error but decible image contains many NaN values.

Code link: https://code.earthengine.google.com/0803e3f7e3518013e310f7d1c666acf9

var fall = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterBounds(table)
.filterDate('2018-10-01','2019-01-01')
.filter(ee.Filter.listContains('transmitterReceiverPolarisation','VV'))
.filter(ee.Filter.equals('instrumentMode','IW'))
.select('VV')
.median()
.clip(table)
.focal_median(50,'square','meters')
.log10().multiply(10.0)
.rename('fall');

Map.addLayer(fall);


texts: Calculating Statistics Based on Expression in QGIS?

texts:
Is it possible to calculate the sum of polygon areas based on whether or not it has a value in another column?

For example:

In column 2 I have the area of polygons. I want to calculate the area of all polygons that contain 'AG' in column 1 or all polygons that are IN column 1. But calculate them separately so I can get the total area of each type in column 1.


texts: Chart number of images within AOI aggregated by month in GEE

texts:
My goal is to chart the number of Sentinel-2 images within AOI aggregated by month on a histogram. Currently I use ui.Chart.feature.histogram and plot by property 'month' that I extracted from image 'system:time_start'. This works for intervals within 1 year, but if I increase the interval it starts agregating corresponding months for different years.

But the goal is to have an X-axis representing Years and Months like on time series chart (below) and Y-axis represent number of images aggregated by month, so how do I do it?

Here is my current implementation:

// ================================ VARS & PARAMS ==========================================

// data

var s2_coll = ee.ImageCollection('COPERNICUS/S2');

// aoi

var aoi = ee.Geometry.Polygon(
  [[[12.073281620012267, 42.652666751867734],
    [12.073281620012267, 42.60214521751987],
    [12.150185916887267, 42.60214521751987],
    [12.150185916887267, 42.652666751867734]]], null, false);

// dates

var year_start = 2019;
var year_end = null; // DO NOT FUNCTION YET

var month_start = 1;
var month_end = 12;

// styling and display

var chart_options = {
    title: 'S2 imagery number',
    hAxis: {title: 'Month'},
    vAxis: {title: 'Image count'},
    colors: ['orange']
};

// misc

// ================================ FUNCTIONS ==============================================

var metadata_collector = function(img) {
  var month  = img.date().format("M"); // get month
  return img.set('month', ee.Number.parse(month)); //
};

// ================================ IMPLEMENTATION =========================================

// preliminary filtering collection by aoi and dates
var s2_coll_filt = s2_coll
    .filterBounds(aoi)
    .filter(ee.Filter.calendarRange(year_start, year_end, 'year'))
    .filter(ee.Filter.calendarRange(month_start, month_end, 'month'));

var s2_coll_filt = s2_coll_filt.map(metadata_collector);

var s2_fc = ee.FeatureCollection(s2_coll_filt);

var num_scn_total = ee.FeatureCollection(s2_coll_filt).size();

var txt_total_scns = 'Total number of scenes within AOI is: ';
var info_text = ee.String(txt_total_scns).cat(num_scn_total);

// ================================ VISUALIZATION ==========================================

// map
Map.centerObject(aoi);

// console
info_text.evaluate(function(result){
  print(result);
});

var histogram = ui.Chart.feature.histogram({
    features: s2_fc,
    // property: 'system:time_start',
    property: 'month',
    minBucketWidth: 1
}).setOptions(chart_options);

print(histogram);

Code Editor script


texts: Saving PDF files as part of QGIS Project

texts:
I am using QGIS to create a record of remediation works across a number of assets. When an asset is remediated, we will often prepare a condition report PDF document which is linked to within the relevant QGIS feature. The condition reports are PDF documents and the link to them can be seen within each feature/data point.
My concern is that if the file path to the PDF document changes, QGIS will no longer be able to identify the PDF or its location.
Is there a way in which the actual report can be uploaded to the QGIS project file? This way, should a folder link change or be amended, the PDF file attached to the QGIS file will not be effected.


texts: How to update attributes of a GeoDataFrame based on another database of the same geometry

texts:
How to update the attribute of a GeoDataFrame based on another database? I have a master GeoDataFrame with line features. I have another GeoDataFrame (say obs) which is a subset of the master geodatabase. My analysis regularly updates the obs. How can I write the results back to the master GeoDataFrame by matching geometry?

Can I do spatial join using geometry as the key?

GeoPandas documentation talks about spatial join by merging. How can I use this approach for updating my attribute values?


texts:
i am trying to develop a webmapping application using openlayers, geoserver, JSF javascript. I published my postgis layers on geoserver and added them to a google basemap on a xhtml page using a javascript code. Now i want to edit one of the layers i published. on the script i added the folowing code :

function  drawpoint(){
          //nouveau_risque = new OpenLayers.Layer.Vector("nouveau risque");
           map.addLayer(nouveau_risque);
           drawControl = new OpenLayers.Control.DrawFeature(nouveau_risque, OpenLayers.Handler.Point);
          map.addControl(drawControl);
          drawControl.activate();
         // drawControl.events.register('featureadded', drawControl, onFeatureAdded);
          drawControl.events.register('featureadded', drawControl, function(f) {

                // create a WKT reader/parser/writer
                var wkt = new OpenLayers.Format.WKT();

                 // write out the feature's geometry in WKT format
                 var out = wkt.write(f.feature);
                 console.log(out);
               });
            // drawControl.deactivate();

           }

           function endDraw(){
               drawControl.deactivate();
           }

and on the xhtml page i wrote :

            <input type="button" value="add point" onclick="drawpoint()"> </input>
            <input type="button" value="stop editing" onclick="EndDraw()"> </input>

so that i can get the geom value of the added point so that i can add it to the database.

The problem is when i run the page, the add control is displayed and i can add a pointy however i don't get the feature's geometry in WKT.

Than you all for your help


texts: Arcpy returning empty string when run in script, but not when run in ArcGIS Desktop

texts:
A part of my script where I extract the view definition is returning an empty string when I run it, but if I run the same exact code in ArcGIS, I get the correct return.

It is only for the desc.whereClause return, the others such as desc.name, desc.nameString are working correctly.

Anyone know what could be causing this?

mxdpath = r"C:\...\blah.mxd"
mapDoc = arcpy.mapping.MapDocument(mxdpath)
dataFrames = arcpy.mapping.ListDataFrames(mapDoc)
for i in range(0,len(dataFrames)):
    lyrs = arcpy.mapping.ListLayers(mapDoc,"",dataFrames[i])
    for j in range(0,len(lyrs)):
        if lyrs[j].isFeatureLayer:
            desc = arcpy.Describe(lyrs[j])
            print "Name string:  " + desc.nameString
            print "Where Clause:  " + desc.whereClause

This code returns a blank u'' for desc.whereClause, however the same code will produce something like 'objectId = 1' if I run it in the python app inside ArcMap Desktop


texts:
I have a line feature that I offset by altering the symbology from "simple line" to "cartographic line" and then adjusting the offset.  However, this segmented the line in many places. I tried a dissolve thinking it was due to being multipart, but it only made a minor change.

Any idea how to correct this?


texts: Difference between QGIS export and gdalwarp

texts:
Does QGIS export/save as use a different projection system/process than gdalwarp? I am working with MASIE ice coverage data available here. My goal is to automate the process of reprojecting the daily images to EPSG 4326; source images are in EPSG 7239. Performing reprojection tests with gdalwarp provided unsatisfactory results, however when the file is manually saved from within QGIS (and providing a CRS of 4326 - see image below) the results are as expected).
My questions are:

Does QGIS Export / Save as CRS projection perform a different reprojection than gdalwarp
To use gdalwarp to achieve the desired output is there an intermediate step(s) required to reproject a raster from EPSG 7239 to EPSG 4326

QGIS Export / Save as with reprojection

QGIS Export / Save as with reprojection results - desired output

gdalwarp code used to reproject
gdalwarp -t_srs EPSG:4326 -r near -of GTiff "/Users/ryangarnett/OneDrive - Global Spatial Technology Solutions Inc/Documents/Data/masie_ice/source/masie_ice_r00_v01_2021005_4km.tif" "/Users/ryangarnett/OneDrive - Global Spatial Technology Solutions Inc/Documents/Data/masie_ice/temp/masie_ice_r00_v01_2021005_4km_4326.tif"

gdalwarp reprojection results - unsuccessful output

UPDATE
Output from gdalinfo based on comment request
Driver: GTiff/GeoTIFF
Files: /Users/ryangarnett/OneDrive - Global Spatial Technology Solutions Inc/Documents/Data/masie_ice/source/masie_ice_r00_v01_2021005_4km.tif
Size is 6144, 6144
Coordinate System is:
PROJCS["Stereographic_North_Pole",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",-80],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-12288000.000000000000000,12288000.000000000000000)
Pixel Size = (4000.000000000000000,-4000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2021:01:06 07:36:19
  TIFFTAG_DOCUMENTNAME=/share/apps/G02186/production/output/geotiff/4km/ice_only/2021/masie_ice_r00_v01_2021005_4km.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.3, Exelis Visual Information Solutions, Inc.
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-12288000.000,12288000.000) (145d 0' 0.00"E, 21d29'36.71"S)
Lower Left  (-12288000.000,-12288000.000) (125d 0' 0.00"W, 21d29'36.71"S)
Upper Right (12288000.000,12288000.000) ( 55d 0' 0.00"E, 21d29'36.71"S)
Lower Right (12288000.000,-12288000.000) ( 35d 0' 0.00"W, 21d29'36.71"S)
Center      (   0.0000000,   0.0000000) ( 80d 0' 0.00"W, 90d 0' 0.00"N)


texts:
I need to develop a QGIS plugin for a project.

First, I launched the following command for the installation on Ubuntu 12.04:

sudo apt-get update
sudo apt-get install qgis
sudo apt-get install python-qgis

Then, I installed the Plugin Builder plugin.

Finally, I tried to follow this simple tutorial.

However, I'm stuck at the execution of the make command, which give to me the following error:

make: pyuic4: command not found

I did some research, but I'm unable to find a solution.
Can someone help me?


texts:
I have installed osm2pgrouting on my Ubuntu 12.04 server with the following code from this tutorial:

sudo add-apt-repository ppa:georepublic/pgrouting-unstable
sudo apt-get update
sudo apt-get install osm2pgrouting

Which ran happily. Unfortunately however, when I go to use osm2pgrouting at the command line, the tool is missing. I realise I could follow the instructions on github but I was hoping to get an easy sudo apt-get install working out of the box.

Thus, how can I get this tool working at the command line? I suspect it is something to do with paths...


texts: ogr2ogr the result of merging two shapefiles has attributes missing

texts:
I am trying to merge two shapefiles from a zip obtained here:
http://wsgw.mass.gov/data/gispub/shape/l3parcels/L3_SHP_M351_YARMOUTH.zip

The two shapefiles contain different attributes for the same town in Massachusetts. I want to join the two shapefiles into one consolidated shapefile with all the attribute information from both files.

I can achieve this on QGIS, however something goes wrong when I try to do this on the command line.

On the command when I do this:

$ ls *.shp
    M351OthLeg.shp  M351TaxPar.shp

$ ogr2ogr -f "ESRI Shapefile" consolidated.shp  M351TaxPar.shp
$ ogr2ogr -f "ESRI Shapefile" -update -append consolidated.shp M351OthLeg.shp

It does produce a consolidated file that is larger than the first two:

$ ls -al *.shp
-rwxr-xr-x@ 1 ssikdar  600    249940 Feb 15  2012 M351OthLeg.shp
-rwxr-xr-x@ 1 ssikdar  600  10765676 Feb 15  2012 M351TaxPar.shp
-rw-r--r--  1 ssikdar  600  11015516 Apr 20 21:55 consolidated.shp

But the new file is missing attributes of the second file.
When I do an ogrinfo:

$ ogrinfo -so -al consolidated.shp
INFO: Open of `consolidated.shp'
      using driver `ESRI Shapefile' successful.

...
SHAPE_Leng: Real (19.11)
SHAPE_Area: Real (19.11)
MAP_PAR_ID: String (26.0)
LOC_ID: String (18.0)
POLY_TYPE: String (15.0)
MAP_NO: String (4.0)
SOURCE: String (15.0)
PLAN_ID: String (40.0)
LAST_EDIT: Integer (9.0)
BND_CHK: String (2.0)
NO_MATCH: String (1.0)
TOWN_ID: Integer (4.0)

It doesn't seem to have the attributes of the last file I tried to merge with it.

$ ogrinfo -so -al M351OthLeg.shp
...
MAP_PAR_ID: String (26.0)
LEGAL_TYPE: String (15.0)
LS_BOOK: String (16.0)
LS_PAGE: String (14.0)
REG_ID: String (15.0)
SHAPE_Leng: Real (19.11)
SHAPE_Area: Real (19.11)
TOWN_ID: Integer (4.0)
TAXPAR_ID: String (18.0)

As a test I merged the two shapefile in Quantum GIS. When I do an ogrinfo on that file, the additional attributes are there:

$ ogrinfo -so -al ~/Documents/m351_qgis.shp
...
MAP_PAR_ID: String (254.0)
LEGAL_TYPE: String (254.0)
LS_BOOK: String (254.0)
LS_PAGE: String (254.0)
REG_ID: String (254.0)
SHAPE_Leng: Real (21.5)
SHAPE_Area: Real (21.5)
TOWN_ID: Integer (10.0)
TAXPAR_ID: String (254.0)
LOC_ID: String (254.0)
POLY_TYPE: String (254.0)
MAP_NO: String (254.0)
SOURCE: String (254.0)
PLAN_ID: String (254.0)
LAST_EDIT: Integer (10.0)
BND_CHK: String (254.0)
NO_MATCH: String (254.0)

Is there something I'm missing in my command line arguments that's giving me a different result than QGIS?


texts: Defining coordinates cells as "Real" using a .csvt file renders them "NULL"

texts:
My .csv file was only adding as a 'legend file' in QGIS, and using the 'Properties' tab I saw QGIS regarding all the entries as text.
In order to resolve this, I made the following .csvt file.

"Integer","String","Real","Real","Integer","Integer"

However, the .csv file still only adds as a 'legend', and using the 'Properties' tab I saw QGIS now regarding the coordinates cells as 'NULL'. All other cells were treated correctly.

How is this possible and how can I resolve this issue?

Sample of my .csv file

id_tourist;name;latitude;longitude;drift;top_list
1;Dam_Square;52,373056;4,892778;1;1
2;Red_Light_District_&_Oude_Kerk;52,3741183;4,8977614;1;1
3;Rijksmuseum_&_Museum_Square;52,3702157;4,8951679;1;1
4;Westerkerk_&_Anne_Frank_House;52,3750066;4,8837106;1;1


texts: Creating "Empty" version of ArcSDE geodatabase?

texts:
Is it possible to create an "Empty" version of an ArcSDE geodatabase?

I have a large default ArcSDE geodatabase that I would like to create an empty version for a new project.  By "Empty" version I am referring to a version that does not contain any of the information maintained in the default version.  As data populated to this "Empty" version it will be post to the default but never reconciled.  The default serves as an overview of all the data I have worked on but new projects do not have any dependencies upon this data (Typically in vastly different geographic locations).

I really would like to utilize only one geodatabase but am open to input and suggestions.

Any ideas?


texts:
I'm trying to loop a ReduceRegion function in Google Earth Engine to obtain a month  mean  NDVI  over a polygon. But I've not been successful yet. I did for one month. At the end, I should have a list or dictionary with 132 mean values of NDVI and date associate.

// Gerar tabela com os valores médios de NDVI mensal da série histórica ????????
// Compute mean NDVI for "limite" for all years. I got it for one month. But don't know to do it for all month in "byMonthYearDoce".

var meanDictionary = ee.Image(byMonthYearDoce.first()).reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: limite,
  scale: 250,
  maxPixels: 1e9
});

print(meanDictionary)

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


texts: sf :: st_intersection() returning duplicate features

texts:
Why does st_intersection add duplicate lines to my sf data frame?

 sf_points

Simple feature collection with 156421 features and 9 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 596640.7 ymin: 7317066 xmax: 607113.9 ymax: 7326187
epsg (SRID):    32722
proj4string:    +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs

sf_polygons

Simple feature collection with 225 features and 4 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 543943.6 ymin: 7317069 xmax: 607110.9 ymax: 7342239
epsg (SRID):    32722
proj4string:    +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs

nrow(sf_points[!duplicated(sf_points), ]) == nrow(sf_points)
> TRUE

st_crs(sf_points) == st_crs(parcelas) # test crs
> TRUE

sf_points_intersection <- st_intersection(sf_points, sf_polygon)

sf_points_intersection
Simple feature collection with 121716 features and 9 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 596646.7 ymin: 7317070 xmax: 607109.9 ymax: 7326184
epsg (SRID):    32722
proj4string:    +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs

nrow(sf_points_intersection[!duplicated(sf_points_intersection), ]) == nrow(sf_points_intersection)
> FALSE

nrow(sf_points_intersection[duplicated(as.data.frame(sf_points_intersection)), ])
> 37056


texts: Getting the status of Google Earth Engine task

texts:
I have queued earth engine task via a Google Service account. The following will generate a dataset and export to Google Cloud Storage
task = ee.batch.Export.table.toCloudStorage(...)
task.start()

I can use task.status() to check the status and the ID of task. Is there any way to get the status of a task for a given task ID?


texts: Dynamic filter on tile (WMS) layer in OpenLayers

texts:
I have a tile layer with a WMS source coming from GeoServer. I want to add and change a filter to it when a user changes their selection. I can see some examples for OpenLayers 3 and Vector layers but the OpenLayers.Filter class has changed and I can't get it to have any effect.

import OlTileLayer from 'ol/layer/Tile';
import OlOr from 'ol/format/filter/Or';
import OlEq from 'ol/format/filter/EqualTo';

this.structureLayer = new OlTileLayer({
  source: new TileWMS({
    url: 'http://test.com/fpat/wms',
    params: {'LAYERS': 'fpat:fpat_surveys',
              'TILED': true},
    serverType: 'geoserver',
    transition: 0
  })
});

 this.structureLayer.filter = new OlOr(new OlEq('structure_types', 'Natural') , new OlEq('structure_types', 'Other'));

When the user changes their selection I'll need to rebuild the filter and make sure the layer updates. The documentation for a vector layer with a spatial filter suggests this, but I I get an error because refresh doesn't exist

this.structureLayer.refresh({force: true})


texts:
In QGIS 3.6, I'm successfully using the "eval (@atlas_pagename)" option that automates the generation of class maps when data
are organized by different fields and making use of the Atlas of the composer option.
However, when trying to use it in the 3D environment the command, it was not recognized.
Am I doing something wrong?


texts: ArcGIS API for JavaScript - Search widget with external geocode service

texts: To be brief, I want to add to my Search widget new geocode service that is not exposed by the ArcGIS Server REST API (more precisely - OneMap REST API - https://docs.onemap.sg/). I wonder if it is possible to add this service as widget source using Locator object and how to do that.


texts: Name points inside a feature made with Geometry.MultiPoint

texts:
i am trying to name the points inside a feature build with a ee. Geometry.MultiPoint function as shown in the image. This points are alredy named in the imports place. But to build this feature a used just lat and long information.

I can´t use a Feature Collection, because this function below needs a feature or a image to run properly.

var featureCollection_getvalue__wet_2000 =
ee.FeatureCollection(getvalue__wet_2000.slice(1).map(function(singleData){
  singleData = ee.List(singleData);
  var dict_getvalue__wet_2000 = ee.Dictionary.fromLists(keys_getvalue__wet_2000, singleData);
  //var point = ee.Geometry.Point([dict.get('longitude'), dict.get('latitude')]);
  //var timeFormat = ee.Date(dict.get('time')).format('YYYY-MM-dd');
  return ee.Feature(geometria, dict_getvalue__wet_2000); //.set('timeFormat', timeFormat);
}));

This function makes it possible to export the results as a list. But without the names in the points the identification of the results is a little difficult.

The results in the exported table are a sequential, but without mentioning who it is at which point.

Does anyone knows who to name the points to see the results with this names?


texts: Create a PostGIS raster using GDAL and Python?

texts:
I'd like to create a PostGIS raster using GDAL and python.

I can get the driver easily:

drv = gdal.GetDriverByName("PostGISRaster")

But how can I create a new raster table? I can't find any information about drv.Create to explain what to put there....?


texts: Displaying multiple quantitative variables without clutter?

texts:
I am working with a set of data which represents total persons at different pay levels. There are 4 different pay grades which ideally would all be displayed on the same map.

I am new to ArcGIS and cartography so don't have many strategies under my belt.

My first strategy was to write a python script to offset each pay grade around the center point and use a label and graduated symbols to express the quantity (It took way too long and looks terribly cluttered)

I also made a map with only one pay grade displayed and it is much clearer although I really want to make a map work with all categories.


texts: How to pad field with blank spaces using calculate field tool of ArcMap?

texts:
I am trying to pad the right side of a field with blank spaces so that when I label the field, it labels further over to the left. The only way that I can seem to make this work so far is by adding a character that arc recognizes, such as a period at the end, otherwise the spaces are not recognized as characters and the label does not move. Any ideas how I can make this work without adding a period or something to the end?

The first image is what it looks like when I label it -  not what I want! The second is what I want to do, label to the left/bottom of the line, but as you can see there is a period at the end of it, which I don't want.

This is the field calculation that I am using, I am sure there must be a better way of doing it.


texts: QGIS GRASS plugin: r.composite quality loss

texts:
When using the r.composite tool within the GRASS plugin in QGIS, the result raster has a cell size of 150x150 while the input rasters have a cellsize of 30x30. Can this be overcome somehow or is this a consequence of the tool?

Data used: Landsat 8, QGIS version: 1.8.0, Operating system: Windows 7.


texts: GeoServer WMS: Retrieved GeoTIFF band contains integer values instead of floats

texts:
I am using gdal_translate to import a GeoTIFF into an oracle database. The import as well as the layer creation in GeoServer works without errors. GeoServer correctly recognizes the band data type as "Real 32 bits".

Unfortunately when I try to retrieve the layer via WMS as GeoTIFF and analyse it the float values are gone and the bands contain only integer values, which is not useful for further image processing.

Any ideas how I can retrieve the bands with their orifinal data types via WMS?


texts:
I have a large dataset of polygons for property sections that I need to reverse geocode. It was recommended that I use a GIS like Quantam to do this however from my research I am not 100% sure if a GIS will support doing this.

The technique that I was going to use until I found GIS was to get the center co-ordinates of the polygon and then pass these co-ordinates into a reverse geo-coding api like google or openstreetmap. However I am not sure if this is very accurate.

The polygon data is in the projection format WGS 84 (EPSG:4326 Lat/Long)

Here is an example of what one of the polygons data looks like:

  MULTIPOLYGON (((174.76394635 -41.34094115,174.763921416699986
  -41.340960566699998,174.763611683300013 -41.340978616699999,174.763597816700013 -41.3409202,174.763793716700008 -41.340691133299998,174.7638963 -41.34071535,174.763958516700001 -41.34073005,174.76394635 -41.34094115)))

What is the best way to reverse geocode this polygon data?

How can I do it on a large scale?


texts:
I want to stack layers of 3 images (median images from 3 seasons). I want to differentiate the land uses in different seasons. for example, there will be more waterbodies in rainy season but less in summer. I want to delineate the disappeared waterbodies in summer and define those as wetlands. The waterbodies remains all through the year will be regarded as waterbody.
I have tried using "add", "subtract" and "mask" function to do that. It didn't work.

My codes:

var maskL8 = function(image_m) {
  var qa = image_m.select('BQA');
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image_m.updateMask(mask);
};

var image_m = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(roi)
.filterDate('2018-06-16', '2018-10-15')
.map(maskL8)
.median());
Map.addLayer(image_m, {bands: ['B5', 'B4', 'B3'], max: 0.3}, 'image_m');

var maskL8 = function(image_s) {
  var qa = image_s.select('BQA');
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image_s.updateMask(mask);
};

var image_s = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(roi)
.filterDate('2018-03-16', '2018-05-15')
.map(maskL8)
.median());
Map.addLayer(image_s, {bands: ['B5', 'B4', 'B3'], max: 0.3}, 'image_s');

var maskL8 = function(image_w) {
  var qa = image_w.select('BQA');
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image_w.updateMask(mask);
};

var image_w = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(roi)
.filterDate('2017-11-16', '2018-02-15')
.map(maskL8)
.median());
Map.addLayer(image_w, {bands: ['B5', 'B4', 'B3'], max: 0.3}, 'image_w');

var stacked = image_m.subtract(image_s);
print('stacked_composite bands', stacked.bandNames());
var stacked_all = stacked.subtract(image_w);
print('stacked_composite bands', stacked_all.bandNames());
Map.addLayer(stacked_all, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'stacked_all');


texts:
I need to estimate the elevation at various access points of some rivers in South Carolina. I use the elevation information to calculate the average gradient (feet/mile) between two access points. The gradient between points A and B on a river is equal to the elevation difference between A and B divided by the distance along the river between A and B.

Usually I get the elevation information from contours on the USGS topo maps. However many South Carolina topo maps have only 20-foot contours when they reach the lower Coastal Plain and within about 20 miles of the coastline. Most states have topo maps with 5-foot contours near the coast.

I learned that South Carolina has been working on acquiring accurate elevation data for the entire state. Digital Elevation Models (DEMs) are available for many counties. The data is available here: http://www.dnr.sc.gov/GIS/lidar.html

They mentioned I would need some kind on software to read the DEM data and  draw the contours with.  They suggested that I try Global Mapper. I tried the trial version but got the following error:
“ Unable to load overlay gdb  -- This release of the GeoDatabase is either invalid or out of date. The ESRI File Geodatabase library does not support GDB files from ArcGIS v9.x  or earlier. You will need to convert your GDB data set to a newer version using ArcGIS v10.x or later to make it readable.”

Next I looked at ESRI’s  ArcGIS Explorer Desktop. I tried to download the trial version, but the larger file was corrupt twice, so I requested a disk be mailed to me.

My question to this group  -- is there another program available to allow me to see the DEM contours? I do not want to pay for an expensive program considering I need relatively little from it.

I would appreciate any advice from users with experience in this area.


texts: How to save raster output in different directory without duplicate in set workspace? Arcpy

texts:
I am using the Weighted Overlay tool in Arcpy and I've not had this issue before with other tools. I'm wondering if it's a tool-specific problem. The output I get from the weighted overlay tool saves in both the workspace I set earlier in the script (this is also where all the layers involved in the overlay are contained) and also the actual output location I want it to go to. I know sometimes arcpy saves temp layers, but it's not saving a temp layer in the set workspace - it does not disappear after closing down the programs. The goal here is to keep the workspace unchanged, and save the output only in the location that I have specified as "outputFolder."

Is there something implicitly wrong in my code?

UPDATE
I have found that the layer being saved in the output folder actually contains data, while the layer being [automatically] saved in the set workspace is empty. From this I'm assuming that the layer in the set workspace is most certainly a "temp" layer but why is it sticking around in my workspace after the process has run and the raster layer as been saved elsewhere?

outputFolder = #some directory here
for stuff in stuffs:
    arcpy.env.workspace = batWrkspace
    #do stuff
    #Create weighted overlay table for Western Washington model
    WstWATable = arcpy.sa.WOTable([[age, spp_weights[0], "VALUE", remap],
                              [dbh, spp_weights[1], "VALUE", remap],
                              [elev, spp_weights[2], "VALUE", remap],
                              [tree_spp, spp_weights[3], "VALUE", remap],
                              [water, spp_weights[4], "VALUE", remap],
                              [bridges, spp_weights[5], "VALUE", remap],
                              [caves, spp_weights[6], "VALUE", remap],
                              [cliffs, spp_weights[7], "VALUE", remap],
                              [mines, spp_weights[8], "VALUE", remap],
                              [talus, spp_weights[9], "VALUE", remap]],
                              [1,3,1])

#Perform weighted overlay for Western Washington model
WO_Result = arcpy.sa.WeightedOverlay(WstWATable)

#Save output
outputName = '{0}_wstwa_hb'.format(folder.lower())
WO_Result.save(os.path.join(outputFolder, outputName))


texts: Is it possible to handle QGis server requests through nodejs

texts:
I have a nodejs based application server. I would also like to add qgis-server functionality to node. Something like this:

expressapp.get(someAuthenticationMiddleware(), qgisServerRequestHandler());

I can't find any reference on how to call/use fast-cgi script through nodejs. Has anyone done it?


texts:
I have two feature classes. Polyline and Point. Both feature classes have the field "NAME" with values populated.

I want to iterate (assuming by field value) with "NAME" field to select polyline and all corresponding points with same "NAME" value that are further than 3 feet from polyline.

I want to create a model in ModelBuilder that will iterate through the entire feature class so i can identify all points that are further than 3 feet away from their corresponding line.

I have successfully added the iterator and added the query in select tool to select both points and polyline by "NAME". However the select by feature tool will not allow the output from select tool to be used as the input feature layer for the select by location tool. It only allows either output to be used for the selection layer.

How do I get the output from the select tool to be used as input feature layer for the select by location tool?


texts: QGIS Layout not exporting although it showed export successfully

texts:
I'm attemptng to export an jpg map from layout using QGIS 3.2.0.
I get the message that the layout has been successfully exported, along with the path, but nothing shows up in the destination specified.
I have tried another computer with QGIS 3.10 (A Coruña) but still not working.
Furthermore, I restarted the program, changed the destination and lower the dpi, nothing helps.

Any ideas?


texts: Clip Raster in ArcGIS Desktop 10.3.1 writes NoData border

texts:
I'm trying to clip a DEM raster using a rectangular polygon shapefile as the extent.

The Clip Raster tool will add/write a NoData border surrounding the output raster. The NoData pixels are displayed as transparent within ArcGIS. When exported as a TIFF for further analysis outside of ArcGIS, there is a black border present.

Does the Clip Raster tool always add a NoData border to any raster output?

I've also tried masking.


texts: OpenLayers or similar map tool with custom travel options?

texts:
I have seen custom fantasy world or game maps using OpenLayers or similar tools.
Now say I want to do exactly that, but I also want to include the in-game traveling options. I imagine it being much like how Google does paths, walk to X, take boat to Y, walk to Z and the like.

Except in games there are other limitations and options. Like instant teleportation to specific points (think fast travel from anywhere), static connection between fixed fasttravel points, one way fasttravel options. And some travel options being locked behind a certain feat in the game (can't fly until you reach X point in the game)
Ideally you'd also take speed of modes of travel into account, so the map can evaluate the quickest path between A and B.

I can't imagine there is a map out there that supports this out of the box or even remotely close to it. But if there were, it would provide support for a wide selection of games. Does anyone know if there is anything like it?


texts: GridMetrics output does not align with base DTM?

texts:
I am trying to derive grid metrics from LiDAR (.las) file using FUSION's GridMetrics tool.

My problem is that the output raster from GridMetrics does not align with the base DTM used to run the tool.

Command:

E:\FUSION\gridmetrics.exe E:\FUSION\testy\N-34-138-B-
b-1-1-4-3.dtm 3 0.5 E:\FUSION\testy\gridmetrics.csv E:\FUSION\testy\*.las

Example:

Point is put on the "offset" pixel center.

Raster with "offset" on coordinates:

Pixel center coordinates are: (628933.5000, 496151.0000)
Pixel left-upper corner coordinates: (628933.2500, 496151.2500)

Raster with "proper" coordinates (and my DTM file that goes into GridMetrics):

Pixel center coordinates are: (628933.7500, 496150.7500)
Pixel left-upper corner coordinates: (628933.5000, 496151.0000)

How can I turn in GridMetrics to output proper pixel coordinates? I also have tried to use switch /extent

E:\FUSION\gridmetrics.exe /extent:E:\FUSION\testy\N-34-138-B-b-1-1-4-3.dtm E:\FUSION\testy\N-34-138-B-
b-1-1-4-3.dtm 3 0.5 E:\FUSION\testy\gridmetrics.csv E:\FUSION\testy\*.las

and /align:

E:\FUSION\gridmetrics.exe /align:E:\FUSION\testy\N-34-138-B-b-1-1-4-3.dtm E:\FUSION\testy\N-34-138-B-
b-1-1-4-3.dtm 3 0.5 E:\FUSION\testy\gridmetrics.csv E:\FUSION\testy\*.las

with my DTM file to no avail.

The raster pixel size is 0.5m. The original DTM was created with LAStools lasgrid.


texts: Writing python script to export map package from each selected record in feature class?

texts:
I am currently building a script to zoom to each providence within a country and export the providence and major cities (I have a layer for each) as a map package.

I need it to select the record, zoom to selected record, run a
definition query, and export the result as a map package with labels as annotations, and repeat for each value within the preselected field.

This is what I have so far:

import arcpy, os

arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

layerName = arcpy.GetParameterAsText(0)
fieldName = arcpy.GetParameterAsText(1)
imgLocation = arcpy.GetParameterAsText(2)
whereList = []

layerFile = arcpy.mapping.Layer(layerName)
layerFile.definitionQuery = ""

print df.name

with arcpy.da.SearchCursor(layerName, fieldName) as cursor:
    for row in cursor:
        for item in row:
            whereList.append(item)

for wl in whereList:
    layerFile.definitionQuery = fieldName + "= '" + wl + "'"
    outFile =  imgLocation + "\\" + wl + ".png"
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToPNG(mxd, outFile, df, resolution=300)

del mxd

Currently the script does export a .png image, however it does not zoom to each selected record. I need some tips for fixing this as well as exporting out each individual result as a map package with annotative labels.


texts:
I´m trying to end this script. It has worked perfectly for other two time range, but for Landsat 5 and Landsat8.

When i tried do the same for LANDSAT7 this error occurred:

  Layer error: reduce.median: 'Can't apply calendarRange filter to objects without a timestamp'.

PART OF MY SCRIPT:

var L7coll = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
  .filterBounds(region)
  .filterDate(startDateL7,endDateL7)
  .map(fDeleteClouds)
  .map(landMask)

//Corrigindo as imagens do Landsat5 a partir de 2003 (que são "riscadas")
  .map(function(image){
                var filled1a = image.focal_mean(2, 'square', 'pixels', 1)
                return filled1a.blend(image);
})

...

var wet_2001 = L7coll.filter(ee.Filter.calendarRange(10,3,'month'));
var wet_2001_median= wet_2001.median();
var dry_2001 = L7coll.filter(ee.Filter.calendarRange(4,9,'month'));
var dry_2001_median= dry_2001.median();

...

Map.addLayer(wet_2001[1],vizParams,'2001_dry_1');
Map.addLayer(dry_2001[1],vizParams,'2001_wet_1');

Any ideas?


texts: GEE: date list to filter an image collection in Google Earth Engine

texts:
I have a CSV with a column of dates. I already imported the CSV, but I would like to know If I could use that column to filter an Image Collection.

Just now, I wrote all dates and used a filter to  exclude dates from the desired list of dates

var Class0Images = image_2.filter(ee.Filter.inList("system:time_start", F0Date));

However, non image was returned and I know there are dates in the list that could return images from the sentinel 1 image collection.

This is the link of my script, which I also computed polarimetric indexes:
https://code.earthengine.google.com/43144b99eb7c14dbdc9bd3953a59dc46

Thanks for your attention and help


texts: How can I add multiple rules in a GetLegendGraphic request in Geoserver?

texts:
I am using following GetLegendGraphcs request

https://develop.blaa.nl/data/test/wms?LAYER=nccw:variabele_woz&request=GetLegendGraphic&format=image/png&styles=nccw_variabele_woz&rule=one

With this I am able to get Legend for Phase 1. I have rules for Phase 2 (rule2) , Phase 3 (rule3) and phase 4 (rule4) too. If I want rules for Phase 1 and Phase 2, how shall I use GetLegendGraphics request where I can include more rules.

Thanks & Regards


texts: Selecting top most feature from overlapping features in ArcGIS Desktop?

texts:
I have an annotation feature class with several overlapping features within the same layer.

While editing in ArcMap 10.3, I am trying to draw a selection box to select multiple features but only want to select the top most features.

Is it possible to limit the selection in this way?


texts: TIN Output empty raster and corrupted Triangulation network

texts:
In Qgis:
I've got a point shape (EPSG:3398) with some height info and want to create a TIN interpolated raster with the Qgis TIN-Interpolation tool.
But after running there was a raster created, which is invisible on the map and with impossible values and the the triangulation network only was made for the southern points like seen in the screenshot.

The point shape doesn't contain corrupted data.

The settings of the TIN tool - no error message while running

The project EPSG is the same than the point shape.
I'm running Q-GIS 3.12.3 on Linux Mint 19.2 (via apt).

What I've tried yet:
Reproject the Point shape to different EPSGs;
Change raster resolution
...with no success.


texts:
I am having trouble with part of my code and it stems from the attributes not updating when I believe they should be updating.

                    if str(f[t1st])==str(f2[start_index]):
                        layer.startEditing()
                        layer.changeAttributeValue(f2.id(),t1st, str(f2[end_index]))
                        layer.changeAttributeValue(f2.id(),t2nd, str(f2[start_index]))
                        layer.commitChanges()
                        print f2[end_index]
                        print f2[t2nd]

This is in the middle of a for loop. The first print statement returns the correct value but the second print value is always NULL. If create a for loop right after, the attributes show as updated and not NULL. Shouldn't the attributes update after I the layer.commitChanges() line? When will the attributes update?


texts: Zonal statistics for millions of polygons

texts:
I have a shapefile with 31 million polygons and I want to extract statistics from a raster for each polygon. Whenever I use zonal statistics I get:

ERROR 999998: Unexpected Error. Failed to execute (Zonal Statistics)

The problem with this error is that there is no clear answer on what is going wrong.

So far I tried breaking the shapefile into smaller pieces and run zonal statistics, changed the processing from background to foreground and changed my working directory to C:/ but I keep getting the same error.

I would like to know what else can I try?


texts: How to access QGIS module from my virtualenv

texts:
I am running an Ubuntu 14.4 and I have installed a python virtual environment.
Then I pip installed several python libraries (gdal, lxml, psycopg2 etc.)

Also in my system I have the module for QGIS.
So when I type python then I can import this module as:

import qgis.core

But when I try to do the same through my virtualenv I get an error: module not found.

When I created the virtualenv I did with the option:

virtualenv --no-site-packages env

Which restricts my env from the system's env. The question is how to access QGIS module in my virtualenv?


texts: Displaying media (JPG) in Popup for Feature Layer?

texts:
var featureLayer = new FeatureLayer({
      url: "https://services.arcgis.com/8MqUfIH3ilrffIf0/arcgis/rest/services/aneityum_SPEED_example_30May2017_geo/FeatureServer",
      outFields: ["*"],
      popupTemplate: {
        title: "{ipa}",
        content: "
Definition: {gloss}
" + "
Definition provided by: {authority}
",
        type: "media",
        mediaInfos: [{
          title: "test",
          //caption: "DEF: {gloss}\n Definition provided by: {authority}",
          type: "image",
          value: {
            sourceUrl: "{image_link}"
            //linkUrl: "{image_link}"
          }
        }]
      }
    });

I'm trying to put images into my popup (using ArcGIS Javascript API), the links to which are included within the GeoJSON file in the map server. However, the media feature does not seem to be working––the "test" title is not even showing up. Any ideas as to why?


texts: Further filtering overpass results set by applying geocodearea has no effect

texts:
I'm trying to apply an additional filter to the results set from an overpass query by simply uncommenting one line, but when I uncomment that line it has no effect and I don't understand why or how to fix it.
Here's the code:
[out:json][timeout:60][bbox:{{bbox}}];

// uncomment this line to further filter
//{{geocodeArea:Augsburg}}->._;

node["tourism"="camp_site"];
out;

way["tourism"="camp_site"];
out center;

relation["tourism"="camp_site"];
out center;

It's line 4 (the geocodeArea line) that I want to uncomment or leave commented out. Here's a link to the script on overpass. Uncommenting line 4 and rerunning the query does not restrict the results.
I realise that I can change line 4 to save the search area into a named result set with:
{{geocodeArea:Augsburg}}->.nameOfMySearchArea;

and that does work but has disadvantages:

I would need to specify the search area's name on the end of all my node, way and relation lines and I'm trying to keep the changes to flip this filter in and out to an absolute minimum (this is important for the intended users).

If I were then to comment out line 4 again, in order to not further filter the results, then it breaks the script and I can find no way to define the search area to something that would not further filter the results. E.g. this does not work:
({{bbox}})->.nameOfMySearchArea;
(Nor a whole load of other things I googled and tried.)

Can anybody advise why the commenting/uncommenting has no effect and what the fix is?


texts: add different custom SVG markers to categorized section using pyqgis

texts: I am trying to find a way to add different custom SVG markers to the categorized symboling. I would like to have different symbols for different values but unfortunately I cannot find a way to do it using script, that's just straight forward through UI.


texts: Aiming to retrieve source and destination from Android, use pgrouting dijkstra to find shortest path and display it back on Android

texts:
I have created an Android application that will accept source and destination from the user. It then calls a php page which takes in the lattitude and longitude of the destination and the source. I am using postgis in the backend, and have created the topology required for pgrouting.
Following is my php code:

First I find the all the source integers thats near to my current location using the following function:-

    function         getSourceId($currentPositionLongitude,$currentPositionLattitude,$approximation,$dbConnectionString){
    // the query will return the source id of the user.
   $db_connection = pg_connect($dbConnectionString);
 $queryForSource = "select source  from mumbai_roads
where ST_Dwithin(geom,ST_GeomFromText('POINT( $currentPositionLongitude $currentPositionLattitude)'),$approximation)";
      $source=array();
        $myresult = pg_query($db_connection,$queryForSource);
        while($myrow = pg_fetch_assoc($myresult)) {
             $resultSource = $myrow['source'];

      array_push($source, $resultSource);

        }
    //print_r($source);
    return $source;
}

Now i find all the target integers that are near my destination location:-

    function getTarget($destinationLongitude,$destinationLattitude,$approximation,$dbConnectionString){
   $db_connection = pg_connect($dbConnectionString);
  $queryForDestination = "select target  from mumbai_roads
where ST_Dwithin(geom,ST_GeomFromText('POINT( $destinationLongitude $destinationLattitude)'),$approximation)";
  $target = array();
    $myresult = pg_query($db_connection,$queryForDestination);
    while($myrow = pg_fetch_assoc($myresult)) {
       $resultTarget = $myrow['target'];

      array_push($target, $resultTarget);
    }
    //print_r($target);
    return $target;
}

Finally as i have both the source and target integers, I try to get the shortest path for all the sources and the targets:-

if(!(empty($source) || empty($target))){
      /*looping through all the source and the targets */
      foreach ($source as $sourceId) {
        foreach ($target as $targetId) {
          $json = getShortestPath($sourceId,$targetId,$dbConnectionString);
          if($json!=""){
            $result=$result.$json;
              echo "$result";
          }
        }

      }
    }

function getShortestPath($source,$target,$dbConnectionString){
   $db_connection = pg_connect($dbConnectionString);
  $queryForShortestDistance = "select seq, id1 AS node, id2 AS edge, di.cost AS cost1, ST_AsGeoJSON(geom) AS geom
  FROM pgr_dijkstra(
    'SELECT id, source, target, cost FROM mumbai_roads',
    $source, $target, false, false
  ) as di
  JOIN mumbai_roads pt
  ON di.id2 = pt.id";
  $geom="";
     $myresult = pg_query($db_connection,$queryForShortestDistance);
    while($myrow = pg_fetch_assoc($myresult)) {
       $geom = $myrow['geom'];

    }
    return $geom;
}

I am getting the path as a multiline which i am mapping on android, but using the above code, i am not getting the whole path, but just a small fraction of the whole path. Please help!


texts: Will a basic pgRouting/OSM data setup allow me to calculate travel-times?

texts: I have a PostGIS database and I need to figure out a quick (in terms of setup) way of calculating the travel-time between points corresponding to geocoded addresses taking into account the route and speed-limits involved. It looks like the pgRouting extension combined with OSM data is going to be my best bet (although please correct me if there's a better alternative). Before I start going down that route, however, I was wondering if someone could just confirm whether or not I'm going to need any additional steps/data/software in order to do what I need? The points I have are within the US states of PA, NJ and DE.


texts: I am looking for the st_dwithin postgis source code. I want to learn the way postgis determines whether a given geometry is within a given radius to another geometry.


texts: GEE linear regression ee.Reducer.linearFit() 'scale'

texts:
I am trying to use ee.Reducer.linearFit() similar to as used in https://developers.google.com/earth-engine/reducers_regression between 'evi' calculated by me and time.

I observed that  absolute values of 'scale' are very high (higher than 1). Since scale corresponds to te slope of the line, I expected the value of it to be in between -1 and +1.

What does high values of scale indicate? And why aren't the scale values a double or float?

Example code here
https://code.earthengine.google.com/d5b329de1a86fed97ca9f53ca71d087c


texts: Moving SQLite database without getting 'bad layers' dialogue in QGIS?

texts:
I would like to move an SQLite db to a new directory on my computer.  The SQLite file is in the same folder as a Qgis project file.  I copied the parent folder to a new location.  When I try to open the Qgis project file, I get a bad layers error for the SQLite db.

What am I doing wrong?

Alternatively, if I try to open the SQLite db from its new location in a new QGIS project, I cannot open it properly.


texts: Is there equivalent to ArcGIS Feature to Polygon tool available?

texts:
Is there a method, combination of tools or a tool which can achieve the same as the Feature to Polygon tool?

I want to get the same result just like I would use that tool.


texts: OSGeo4w installer fails: unable to extract /appa/qt5/bin/qt5webenginecored.pdb

texts:
I am trying to install OSGeo4w. I have plenty of RAM and disk space, running Server 2016 standard. Installation fails with the following log:

2019/07/08 09:13:30 Menu name: OSGeo4W
2019/07/08 09:13:30 net: Direct
2019/07/08 09:13:34 site: http://www.norbit.de/osgeo4w/
2019/07/08 09:15:50 Extracting from file://C:\Temp/http%3a%2f%2fwww.norbit.de%2fosgeo4w%2f/x86_64/release/qt5/qt5-libs-pdb-debug/qt5-libs-pdb-debug-5.9.2-1.tar.bz2
2019/07/08 09:17:48 Failed to output cygfile:///apps/Qt5/bin/Qt5WebEngineCored.pdb
2019/07/08 09:20:19 note: Installation Complete
2019/07/08 09:20:19 Ending OSGeo4W install


texts: How to set a spatial filter with Python/GDAL?

texts:
Is there a way with GDAL (Python) to set a spatial filter (bounding-box) on a raster before I process it? I want to read it as an array.

BandReadAsArray(band, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None, buf_ysize=None, buf_obj=None)


texts: OpenLayers Control.WMSGetFeatureInfo does not send request to geoserver

texts:
I'm refactoring the code used to getfeatureinfo from geoserver running on port 8080 at http://www.edpp.cz/zul_mapa-povodnoveho-planu-mesta/. The old version uses OpenLayers.loadURL and successfully requests the response from the server. The code is like:

map.events.register('click', map, function (e) {
        var params = {
                REQUEST: "GetFeatureInfo",
                EXCEPTIONS: "application/vnd.ogc.se_xml",
                BBOX: map.getExtent().toBBOX(),
                SERVICE: "WMS",
                VERSION: "1.1.1",
                X: e.xy.x,
                Y: e.xy.y,
                INFO_FORMAT: 'text/html',
                QUERY_LAYERS: layers_on,
                FEATURE_COUNT: 1,
                Layers: layers_on,
                Styles: '',
                WIDTH: map.size.w,
                HEIGHT: map.size.h,
                format: format,
                srs: map.layers[0].params.SRS};

            if(offline)
                OpenLayers.loadURL("http://localhost:8080/geoserver/wms", params, this, setHTML, setHTML);
            else
                OpenLayers.loadURL("http://www.edpp.cz/mapinfo/geoserver/wms", params, this, setHTML, setHTML);

        OpenLayers.Event.stop(e);
    });

With the new code I try something like this:

popup = new OpenLayers.Control.WMSGetFeatureInfo({
            url: offline ? 'http://localhost:8080/geoserver/wms' : 'http://www.edpp.cz/mapinfo/geoserver/wms',
            layers: unbaseLayers,
            layerUrls: offline ? ['http://localhost:8080/geoserver/wms'] : ['http://www.edpp.cz/mapinfo/geoserver/wms'],
            title: 'Identify features by clicking',
            queryVisible: true,
            maxFeatures: 1,
            eventListeners: {
                getfeatureinfo: function(event) {
                    if (console) {
                        console.log(event);
                    }
                    if (event.text.length > 1088) { // nezobrazí prázdné bubliny
                        map.addPopup(new OpenLayers.Popup.FramedCloud(
                            "popup",
                            map.getLonLatFromPixel(event.xy),
                            null,
                            event.text,
                            null,
                            true
                        ));
                    }
                },
                nogetfeatureinfo: function(event) {
                    if (console) {
                        console.log(event);
                    }
                }
            }
        });

If I look into Firebug network log, I see that the clicking the map doesn't send any request to the geoserver (it runs fine on localhost though). Is this just a proxy related problem? I've already read about need to set up a proxy when running geoserver on the different port than a website. But it's really puzzling me there is no request being sent. Might it be a typo? I can't see any.

What if I am not able to set up a proxy, do I have any other option?

PS: All I get as a result is object of nogetfeatureinfo type.


texts:
I'm a newbie to POSTGIS and I wonder if someone could help me out.

I have loaded all postal codes with their geographical coordinates (geography) and created an index on this coordinate.

When I run the following query it works just fine (within a few millisecond over a few million instances)

SELECT
  postcode, huisnummer, huisletter, huisnummertoevoeging,
  ST_X(geolocatie :: geometry) as latitude,
  ST_Y(geolocatie :: geometry) as longitude,
  st_distance(geolocatie, ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326)) as afstand
FROM inspire
WHERE ST_DWithin(geolocatie, ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326), 1000)
ORDER BY afstand
LIMIT 5;

geolocatie is of type geography and I have defined the following index on geolocatie

CREATE INDEX inspire_geolocatie_index
    ON public.inspire USING gist
    (geolocatie)
    TABLESPACE pg_default;

This works fine as long as there are neighbors, but the query seems to have difficulties for the use case that there are no neighbors. Any suggestions please?


texts: Run module function over FeatureCollection and write results into table in GEE

texts:
This is a follow up question on (Extract and plot from an Array in GEE) but can also be treated as a standalone question.

I want to run a LandTrendr analysis for a given number of points (coordinates) and write the results into a table such as

feature.id | year | original | fitted
00_00      | 2001 | 1001.00  | 1002.10
00_01      | 2002 | 1021.00  | 1050.20
...
01_00      | 2001 | 2001.00  | 2002.10
01_01      | 2002 | 3001.00  | 3002.10
...

My function over the FeatureCollection is not working and I can't get it right.

// Points
var p1 = ee.Geometry.Point([-122.64550158728298,38.44643840127136]);
var p2 = ee.Geometry.Point([-123.06960473413164,38.54089331487732]);
var p3 = ee.Geometry.Point([-122.85149523962673,38.33988067525676]);
var fg_points = ee.FeatureCollection(ee.List([ee.Feature(p1),ee.Feature(p2),ee.Feature(p3)]));

Map.setCenter(-122.64550158728298,38.44643840127136, 11);
Map.addLayer(feat);

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

// define parameters
var startYear = 1985;
var endYear = 2018;
var startDay = '06-20';
var endDay = '09-20';
var index = 'NBR';
var ftvList = ["NBR"];
var maskThese = ['cloud', 'shadow', 'snow'];

//// function to map over the FeatureCollection
var mapfunc = function(feat) {

  var geom = feat.geometry()

  var collection = ltgee.buildSRcollection(startYear, endYear, startDay, endDay, geom, maskThese)
  var timeSeries = ltgee.buildLTcollection(collection, index, ftvList)

  var lt = ee.Algorithms.TemporalSegmentation.LandTrendr({
    timeSeries:             timeSeries,
    maxSegments:            6,
    spikeThreshold:         0.9,
    vertexCountOvershoot:   3,
    preventOneYearRecovery: true,
    recoveryThreshold:      0.25,
    pvalThreshold:          0.05,
    bestModelProportion:    0.75,
    minObservationsNeeded:  6
  });

  var fittedImage = ltgee.getFittedData(lt, startYear, endYear, index);

  var fittedCollection = ee.ImageCollection(fittedImage.bandNames()
  .map(function (yearString) {
    yearString = ee.String(yearString);
    return fittedImage.select(yearString)
      .set('system:time_start', ee.Date.parse('yyyy', yearString).millis())
      .rename('Fitted');
  })
  );

  // get feature id
  var id = ee.String(feat.id());

  // make an empty list to store the features
  var newfc = ee.List([])
  //// function to iterate over fittedCollection
  var addVals = function(img, fc) {
    // the initial value is the empty list
    fc = ee.List(fc)
    // get the date as string
    var date = img.date().format()
    //// extract the value of 'Fitted'
    var value = img.reduceRegion(ee.Reducer.first(), geom, 30).get('Fitted')
    // If the value is null then store it as 'No data'
    var val = ee.String(ee.Algorithms.If(value, ee.String(value), ee.String('No data')))
    // make the name of the feature (feat_id-date)
    var featname = ee.String("feat_").cat(id).cat(ee.String("-")).cat(date)

    // make the Feature
    var newfeat = ee.Feature(geom, {name:featname,
                                    value:val})
    // add the value to the list
     return fc.add(newfeat)
  }
  var newfeat = ee.FeatureCollection(ee.List(yearly.iterate(addVals, newfc)))
  return newfeat;
};

var newft = fg_points.map(mapfunc).flatten();

Export.table.toDrive(newft,
"export_Points");

https://code.earthengine.google.com/1acf553943e0bba49bca66878e8644c8


texts: Converting WGS to latitude and longitude in QGIS?

texts: How do you convert the WGS coordinate to latitude and longitude?


texts: Tilemill/CartoCSS: how to center label in a polygon

texts:
For some reason, labels are always placed in the bottom right of a polygon. I want them in the center of a polygon.

I tried these:

text-horizontal-alignment: middle;
text-vertical-alignment: middle;

but it doesn't help.

Update:

I noticed a strange thing: buildin labels are centered inside a polygon, but for regular area=yes or landuse=yes they are along the line, i.e. text-placement:interior is not working consistently.

Full example:

#labels {
    text-name: "[name]";
    text-fill: #444;
    text-face-name: 'Alfios Bold';
    text-placement: interior;
    text-size: 12pt;
}

#cemetery[building='yes'] {
    polygon-fill: #0f0;
}

#cemetery[landuse='residential'] {
    polygon-fill: #0f0;
}


texts:
Is there a way to preserve the format of the legend in an .lyr file in ArcMap?

The legend is simple - bold layer name atop the symbol with no label or description. When I save it as a .lyr file, it instead imports with the header and label and does not keep the layer name at the top.


texts: Can I upload a shapefile to the GeoServer Import Data page, or does it only import data locally on the server?  When I click Browse to select the data to import, it seems to only browse the server filesystem.


texts: Vector Fields as Input Variable in Field Calculator

texts:
I’m working with QGIS 3.2 Bonn and I try to build a workflow in the graphic modeler. The model aims to build variable buffers depending on values from two attribute fields of an input vector layer. For calculating the distances for the dynamic buffers I use the field calculator and the input fields to the field calculator are addressed by their distinct names.

This only works as long as the attribute fields of the input layer always have the same names, so I wanted to modify the model for making it more flexible to other input vector layers with different field names. For this I tried to use the Vector-Field variable as input to the field calculator. In the Variable distance buffer tool, I managed to use one Vector-Field as input variable, but with the field calculator I haven't been succesfull yet.

Is it possible to use Vector Field as input in the modeler, or do I have to write a script for that?

Would be really glad if you could help me.

(I also searched in gis-stackexchange and the net but couldn’t find an answer to this for QGIS3.x.
I found only other questions with similar problems but without any solution:
Graphical modeler in QGIS 3.0.3 - Access Vector Field in Calculator
or questions refering to QGIS2.x:
How to link Number value with Field calculator in QGIS Modeler ?
)

An extract of my model with the part I don't get worked:

Here a screenshot of the field calculator and my idea of how to use it:


texts: Run PyQGIS from outside QGIS on MacOS HighSierra

texts:
I know this question has already been asked multiple times but I just can't figure it out.
I am working on MacOS HighSierra 10.13.6 and QGIS 3.16. I also have Python 3.6 and 3.9 installed. Python 3.6.7 is accessed via python3 in the terminal.
I try to run QGIS algorithms within my python script and am therefore lookin for sth like this:
import os
import sys
sys.path.append('/Applications/QGIS.app/Contents/Resources/python/')
sys.path.append('/Applications/QGIS.app/Contents/Resources/python/plugins')
from qgis.core import *

But I always get this error:
ImportError: dlopen(/Applications/QGIS.app/Contents/Resources/python/qgis/_core.so, 2): Library not loaded: @rpath/qgis_core.framework/Versions/3.16/qgis_core
  Referenced from: /Applications/QGIS.app/Contents/Resources/python/qgis/_core.so
  Reason: image not found


texts: A numerical sequence in Qgis attribute table

texts: In my attribute table of a shapefile I created a new field ('number') and I would like to fill it with: