texts: Keeping only highest point in cluster of points using ArcMap

texts:
i have a cluster of points in particular areas around NZ.
I am only wanting the point with the highest value in this cluster to remain. Also many of the points have the same value, i am just wanting 1 point in the cluster to remain. It does not matter which highest point it is, as long as it is one with the highest value.
Any idea on how to do this?

Bottom pic is using focal stats > raster calc and then raster to point


texts: Can SpatiaLite Indexes be used from views?

texts:
I'm working with a large SpatiaLite database and trying to take advantage of some spatial indexes I have on my tables that I'm attempting to access through views.

Essentially, I have a number of views which each query multiple spatial tables and union the results. Before I tried using the views, I was running queries against the tables themselves and was able to use the spatial indexes in a similar manner as follows:

SELECT lc1.lc_name AS "Local Council",
  lc2.lc_name AS "Neighour"
FROM local_councils AS lc1,
  local_councils AS lc2
WHERE lc2.ROWID IN (
  SELECT pkid
    FROM idx_local_councils_geometry
    WHERE pkid MATCH RTreeIntersects(
      MbrMinX(lc1.geometry),
      MbrMinY(lc1.geometry),
      MbrMaxX(lc1.geometry),
      MbrMaxY(lc1.geometry)));

[From the spatialite cookbook]

As I understand it, this is the only way that spatial indexes can be used (i.e., not implicitly) due to a SQLite limitation. If this is the case, then am I correct that views cannot take advantage of spatial indexes? I would imagine that the bounding box information would have to be hard coded into the queries in the view or entered dynamically somehow when the query against the view is made.

Am I missing anything or are these assumptions correct and my best option is to stick with my individual queries against the tables themselves?


texts: How to define a custom function that returns value of a field in an unsaved form

texts:
Filter expression builder in QGIS 2.8 allows us to define custom functions using python. I'm intrested would it be possible to define a function which returns a value of a field in an unsaved form that is currently being filled. I would like to use this output in a filter expression for value relation.

(see: Limiting available values in dropdown field of QGIS form according to value selected in previous field?)

Possibility of reading the contents of unsaved fields is mentioned in theese two posts, but i lack the knowledge of python to adapt it to purpose mentioned above:

QGIS Value Map Edit Widget with properties dependent on another field

How to automatically populate fields instantly?

If this is not overly complicated i would be very gratefull for an example of code i could addapt to my purpose.


texts: Error message when using v.in.ascii in GRASS

texts:
I'm trying to follow the this GRASS tutorial and I'm already stuck at the first step. I have acii points in an .xyz and supposed to import them with v.in.ascii. However, I get the following error message

  WARNING: Expected 5 columns, found 4 columns
  ERROR: Broken row 384: '32346000.00 5630383.00  100.01'
  (Sun Feb  2 16:02:14 2020) Befehl ausgeführt (0 Sek)
  (Sun Feb  2 16:05:43 2020)
  v.in.ascii -n input=C:\Users\denni\Desktop\xyz_daten\dgm1.xyz output=dgm1 separator=space text=space z=4

I don't find the setting to tell the algorithm to expect only 4 columns.

That's the content of the ascii file:

And that's the settings under the "point" field:

how can I solve the error?


texts: Encoding problem with Portuguese characters while opening shapefile

texts:
I use a shapefile with names in Portuguese that contain special characters. When I add it to QGIS it exchanges special characters for symbols, so I have problems when comparing values with other layers.

I use QGIS - 3.8.3 Zanzibar on Ubuntu 18.04. So how can I correctly encode the special characters in the shapefile?


texts: Connect ESRI Feature Service to Oracle Spatial database

texts:
I need to display some business and feature data from an Oracle Spatial database on a website. The website will use ESRI's HTML5 map components (https://developers.arcgis.com/en/javascript/jssamples/#search/html5), which would get their features from some kind of Feature web service.

So I'm looking for a web service by ESRI which can be hooked on an Oracle Spatial data source. I don't know if I need a WMS, WFS or an OGC service (sorry, no experience with such services), as long as it's compatible with ESRI's HTML5 Javascript libraries.

Is there such a web service component in ESRI's software suite? Or do I need the full ArcGIS Server (including ArcSDE) to display Features in a website?

What kind of web service (WMS, WFS, OGC, etc.) would best be suited for such a scenario?


texts: Can I store raster data in a spatialite database (in QGIS)?

texts:
This might be a super simple question, but I'm a little confused. I don't know what the best way is to share qgis-projects that have raster-data, vector-data and some applied styling. I figured a database would be a good option. And as a spatialite-database can be shared easily, that might be a good option.
But I don't manage to put a raster into it. I think there is some extension-driver called RasterLite, but I also don't know how this works properly...


texts:
I have a series of MapInfo layouts that I use as templates for various standard map sizes (eg A0, A4 etc) that contain title blocks, and I wondered what is the best method for getting these into QGIS, if possible.

If not, can someone suggest how to create standard frames and title blocks?


texts:
I am new in this field and I am having a problem with my DEM, (SRTM 30M).

When I calculate flow accumulation a straight line appears if I fill the DEM. And rivers are not connected if I don't fill it.

How can I solve this problem?


texts: Load layers from SQLite database in Model Designer

texts:
Regularly, I have to load layers from SQLite databases through Spatialite in Qgis 3.16 and perform some algorithms. Extension of these SQLite databases is not *.sqlite but another format (proprietary).
My step 1) is to load layers from SQLite database in Qgis through Spatialite.
My step 2) is the use of models to perform algorithms on these layers.
I would like to load databases layers directly from models algorithms, thus to avoid this manual step of loading layers with Spatialite.
I tried many ways to reach this goal, and finally circumvented this issue. Inspired from Load lines from .OSM in Processing Modeler, I am now able to load temporary layers with the appropriate CRS from SQLite database into Qgis directly. Here is the scheme of the model:
Due to the "String concatenation" algorithm, I may not be able to use this model scheme to load layers for some databases (the names of the layers could slightly change). Using the QGIS documentation  and Stack Exchange Q&As, I am now able to do some SQL Queries with DB Manager. I think that the use of SQL Queries could solve this issue, also may be easier and faster. However, I did not find an algorithm or a way to execute AND load SQL from a SQLite database, as the algorithm "PostgreSQL execute and load SQL" do for Postgis database.
If I understood it correctly, I also need to find a way to declare in an SQL query:

an Unique ID field name

a Geometry field name

for each layer. Is there a way to do it? Please be clear, I'm still not very comfortable with SQL Queries.


texts: Check if a point is inside a polygon from .shp file

texts:
I'm trying to build an application that will resolve the corresponding "region" of a point. So basically if a lat-lng point is inside Polygon-A or Polygon-B.

I have all the files (I hope at least), a .shp, a .prj and two more binary files.
I'm trying to do it in Go, I've found a library that is able to read the shp file (go-shp: https://github.com/jonas-p/go-shp )  but I'm not sure how to derive the coordinates.

Is it possible or am I missing something? Do I need to convert it in GeoJSON or something like that?

Sorry but I'm pretty noob about GIS!

EDIT:
Reading the Polygon object returned by the library the points associated are in the form of something like this:

[
    {984002.7992000002 4.6707144145e+06}
    {985067.0396999996 4.670245166999999e+06}
    {985239.5382000003 4.670172667400001e+06}
    ...
]

and I'm not sure on how to use it.

The .prj file contains these informations:

PROJCS[
    "WGS_1984_UTM_Zone_32N",
    GEOGCS[
        "GCS_WGS_1984",
        DATUM[
            "D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]
        ],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",9.0],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]
]


texts: Password protect network link on Geoserver

texts: I plan on serving data with geoserver to multiple people in my organization. The data is somewhat sensitive and I was wondering if there is a way to password protect the info so that in order to visualize it on Google Earth as a Network link, you would need a global password?


texts: Google Earth Engine with Folium and Sentinel-1 Plotting

texts:
I'm attempting to plot Google Earth Engine Sentinel-1 SAR data in a Jupyter notebook with Python and Folium.  To do this I'm trying to combine two pieces of documentation:

EE Python Install -> Map Visualization -> Interactive Map
Sentinel-1 Algorithms -> Metadata and Filtering (JavaScript)

What I've pieced together is currently freezing at the getMapId(vis_params) stage:

import ee
import folium
ee.Initialize()

# Load the Sentinel-1 ImageCollection
sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD')

               # Filter to get images with VV and VH dual polarization.
vh = (sentinel1.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
               .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
               # Filter to get images collected in interferometric wide swath mode.
               .filter(ee.Filter.eq('instrumentMode', 'IW')))

# Filter to get images from different look angles.
vhAscending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
vhDescending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))

# Create a composite from means at different polarizations and look angles.
composite = ee.Image.cat([
  vhAscending.select('VH').mean(),
  ee.ImageCollection(vhAscending.select('VV').merge(vhDescending.select('VV'))).mean(),
  vhDescending.select('VH').mean()
]).focal_median()

# Display as a composite of polarization and backscattering characteristics.

# Define a method for displaying Earth Engine image tiles to folium map.
def add_ee_layer(self, ee_image_object, vis_params, name):
  map_id_dict = ee.Image(ee_image_object).getMapId(vis_params) # STUCK
  folium.raster_layers.TileLayer(
    tiles = map_id_dict['tile_fetcher'].url_format,
    attr = "Map Data © Google Earth Engine",
    name = name,
    overlay = True,
    control = True
  ).add_to(self)

# Add EE drawing method to folium.
folium.Map.add_ee_layer = add_ee_layer

# Fetch an elevation model.
dem = ee.Image('COPERNICUS/S1_GRD')

# Create a folium map object.
my_map = folium.Map(location=[20, 0], zoom_start=3, height=500)

# Add the elevation model to the map object.
my_map.add_ee_layer(composite, {'min': [-25, -20, -25], 'max': [0, 10, 0]}, 'composite')

# Add a layer control panel to the map.
my_map.add_child(folium.LayerControl())

# Display the map.
display(my_map)

Would anyone who has experience with this be able to point out what I'm doing wrong?


texts: Creating Scale field based on polyline length (ArcMap)(Python)

texts:
I'm playing around with data driven pages, and have a polyline feature class that I'm using as my index. Now, because the lines are all of variable lengths I can't just use the default settings as some pages will be way too zoomed out. Instead I want to create a scale field to use, calculated based on the length of the line. What I have found so far is how to calculate it based on the difference between the  length and height but this doesn't quite give the required results:

import arcpy
from arcpy import env
import os

env.workspace = r"U:"
mxd = arcpy.mapping.MapDocument(r"U:\EMS.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
newWeed = r"U:\ArcGIS\EMS\split_weed_index.shp"

print "Adding SCALE field"
addField =[['SCALE','TEXT', "", "", "", "", 'NULLABLE', "", ""]]
for item in addField:
    arcpy.AddField_management(newWeed, item[0])

scalefield = "SCALE" #scale field added above
lengthfield = "L" #actual length field

print "Calculating scales"
with arcpy.da.UpdateCursor(newWeed ,["SHAPE@", scalefield]) as cursor:
    for row in cursor:
        ext = row[0].extent
        maxdiffxy = max([ext.XMax-ext.XMin,ext.YMax-ext.YMin])#max width or
height of feature
        if maxdiffxy<1000:
            row[1] = 6000
        else:
            row[1] = 10000
        cursor.updateRow(row)

print "Scales calculated"

print "Done"

Ideally I'd want something along the lines of "if the length is 1700m or more set the scale to 1:10000, otherwise 1:6000


texts: Load Raster from MS SQL Server into QGIS

texts:
Attempting to connect to an MS SQL database in QGIS.
Database connects fine, and all the vector items load correctly.
Rasters, however, appear as polygons, geometry column shows as FOOTPRINT, and when loaded in they are a basic polygon of the raster footprint, not the raster file.
It's from an organisational server, set up for use in Arc as an SDE, but personal requirements mean working in QGIS would be better. The rasters are  in SDR format - is it a case this isn't supported in QGIS?
Anyone able to advise on how to load the rasters from the MS SQL Server?

The system is already pre-set up by our company in Arc 10.3.1 - we navigate in the catalog to a location where there are  sde files - Properties in Arc describe it as a SQL Server Enterprise Geodatabase. The connection properties give the server instance
Loaded into QGIS 3.4 via the datasource manager using the server instance as the host, selecting the required database, and have tried testing the 4 database options with no luck


texts: Create NetCDF file with irregular rotated lat-lon grid for QGIS

texts:
I want to write data to a NetCDF file with an irregular rotated lot-lan coordinate grid and open it in QGIS. The data is on an irregular grid which follows a satellite swath. The grid is irregular in the sense that the distances between points can vary. Basically, I have a matrix of wind speed values ws10 at certain latitude and longitude coordinates. How do I properly save these ws10(lat, lon) to a NetCDF file, so that they can be viewed in QGIS projected on a world map?
This is what my NetCDF file header currently looks like:
netcdf waves {
dimensions:
        y = 52 ;
        x = 60 ;
variables:
        float lat(y, x) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        float lon(y, x) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;
        float ws10(y, x) ;
                ws10:long_name = "Wind Speed at 10 m" ;
                ws10:units = "meters per second" ;
                ws10:coordinates = "lon lat" ;
}

If I open the file in QGIS and add it as a raster or mesh layer (based on Dave X's comment) it takes the x-y-indices as coordinates. For some reason the grey-scale raster layer is flipped:

Yet Panoply projects the raster as intended:

Based on other questions, I have tried to add the CRS information to the file as shown below, but QGIS still uses the x-y-grid. How can this be solved?
netcdf test_man {
dimensions:
        y = 52 ;
        x = 60 ;
variables:
        int crs ;
                crs:spatial_ref = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.017453292519943295]]" ;
                crs:grid_mapping_name = "latitude_longitude" ;
                crs:longitude_of_prime_meridian = 0. ;
                crs:semi_major_axis = 6378137. ;
                crs:inverse_flattening = 298.257223563 ;
        float lat(y, x) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        float lon(y, x) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;
        float ws10(y, x) ;
                ws10:long_name = "Wind Speed at 10 m" ;
                ws10:units = "meters per second" ;
                ws10:coordinates = "lon lat" ;
                ws10:grid_mapping = "crs" ;
}


texts: How to convert excel csv file to filled contour Google Earth Map?

texts:
I have a csv file http://www.geocities.ws/radar/data2.csv.
The first column is latitude.
The second column is lonitude.
The third column is precipitation type.
5 represents freezing rain. 15 represents sleet. 25 represents rain. 35 represents snow. 45 represents hail.
How would I translate that into a filled contour plot like the following? (http://www.mathworks.com/help/examples/graphics/CreateAFilledContourPlotExample_01.png)
I would like this on a Google Earth kml file. How would I edit the spacing between the contours?


texts: Why is QGIS (resource sharing plugin) telling me that R (3.4.4-1ubuntu1) is not installed?

texts:
This morning I went to install the R Processing scripts to my installation of QGIS3.12 (on linux Mint 19.3) and although the Resource Sharing plugin installs the resources required and the scripts appear in the Processing pane, I am told that R is not installed when it is (first image below)
The dialog that appears is in the second image.
Relatedly, the third image below shows my Options for the R provision with an altered R directory (as per comments' suggestions). When these settings are confirmed, QGIS crashes and leaves the following crash report:
The package providing /newroot/app/bin/qgis could not be found.


texts: How to extract DEM (Digital Elevation Model) from 3D model

texts: How can I extract a Digital Elevation Model from a 3D model of a building, made with SketchUp, using Blender,qgis or any other software?


texts:
I am trying to update the geometry of a point FC stored in a FileGeodatabase using ArcPy.

I got this Python script to compare the values of the SHAPE@X and SHAPE@Y and change the Geometry if not equal to the values stored in the attribute table (X_UTM_WGS84, Y_UTM_WGS84).

Even if I edit the feature and remove the extra digits past the second digit I get this when I print the values: I do not understand why the values are different ???

***** CODE USED:

2015-02-26_Compare_Geometry_With_Table_XY.py
 Compare the Geometry of points using the table values

import arcpy

fc = "D:/DATA_2014.gdb/UTM/TEST_PROJECTION_PT"
fields = ('SHAPE@X', 'SHAPE@Y', 'X_UTM_WGS84', 'Y_UTM_WGS84', 'POINT_NAME')

 Create search cursor for feature class

i=0
with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:  # row[4]+ POINT NAME OUT OF PRINT
        i+=1
        print 'POINT ---> ' + str(i)
        print '\tSHAPE@X= ' + "%.11f" % round(row[0],11) + ' X_UTM_WGS84 = ' + "%.11f" % round(row[2],11) + ' EQUAL= ' + str(row[0] == row[2]) + ' DIFF= ' + str(row[0] -row[2])
        print '\tSHAPE@Y= ' + "%.11f" % round(row[1],11) + ' Y_UTM_WGS84 = ' + "%.11f" % round(row[3],11) + ' EQUAL= ' + str(row[1] == row[3]) + ' DIFF= ' + str(row[1] -row[3])

del cursor, row
print '\t\tEND....'

***** RESULTS:

POINT ---> 1
    SHAPE@X= 460965.17999999970 X_UTM_WGS84 = 460965.17999999970 EQUAL= True DIFF= 0.0
    SHAPE@Y= 1288230.55000000075 Y_UTM_WGS84 = 1288230.55000000005 EQUAL= False DIFF= 6.98491930962e-10
POINT ---> 2
    SHAPE@X= 413207.26999999955 X_UTM_WGS84 = 413207.26999999955 EQUAL= True DIFF= 0.0
    SHAPE@Y= 1382676.34999999963 Y_UTM_WGS84 = 1382676.34999999963 EQUAL= True DIFF= 0.0
POINT ---> 3
    SHAPE@X= 433224.66000000015 X_UTM_WGS84 = 433224.65999999997 EQUAL= False DIFF= 1.7462298274e-10
    SHAPE@Y= 1363812.55000000075 Y_UTM_WGS84 = 1363812.55000000005 EQUAL= False DIFF= 6.98491930962e-10
POINT ---> 4
    SHAPE@X= 428742.87000000011 X_UTM_WGS84 = 428742.87000000000 EQUAL= False DIFF= 1.16415321827e-10
    SHAPE@Y= 1233218.47000000067 Y_UTM_WGS84 = 1233218.46999999997 EQUAL= False DIFF= 6.98491930962e-10
POINT ---> 5
    SHAPE@X= 407040.38999999966 X_UTM_WGS84 = 407040.39000000001 EQUAL= False DIFF= -3.49245965481e-10
    SHAPE@Y= 1227264.44999999925 Y_UTM_WGS84 = 1227264.44999999995 EQUAL= False DIFF= -6.98491930962e-10
POINT ---> 6
    SHAPE@X= 403344.58999999985 X_UTM_WGS84 = 403344.59000000003 EQUAL= False DIFF= -1.7462298274e-10
    SHAPE@Y= 1250201.60999999940 Y_UTM_WGS84 = 1250201.61000000010 EQUAL= False DIFF= -6.98491930962e-10
POINT ---> 7
    SHAPE@X= 420615.79999999981 X_UTM_WGS84 = 420615.79999999981 EQUAL= True DIFF= 0.0
    SHAPE@Y= 1249885.83999999985 Y_UTM_WGS84 = 1249885.83999999985 EQUAL= True DIFF= 0.0
        END....


texts: Seeking common interface for ST_Geometry SDO_Geomtry

texts:
I've worked a fair bit with Oracle Spatial and I've had some requests to transition to PostgreSQL Spatial.

Is there a common interface for SDO_Geometry and ST_Geometry?  I'm working in Java with EclipseLink ORM.


texts: OpenLayers WFS-T tinyows insert transaction: error GID value

texts:
I try to implement insertion of new point using Openlayers v2.13 saveStrategy calling tinyows which implements the WFS-T transaction; I use the example from http://mapserver.org/tinyows/openlayershowto.html  - the code is exactly the same except for the layer.

I have two questions:

1) adding a new point implements a transaction which uses "MULTIPOINT" but my postgis layer is of type POINT - is there a way to force "POINT" ? My solution up to now is to use a TRIGGER in postgresql to do an INSTEAD OF to change MULTIPOINT TO POINT.

2) I get an error because the insertion creates a "gid" value which should not be included, as it is added by default by the counter. The TINYOWS log give me the following transaction and error:

SQL

  INSERT INTO "arch_geo_dbt"."p030104_v" ("gid","the_geom") VALUES ('042159162153000033','010400002031BF0D00010000000101000000B69CF08B3069314130B8664669AD5541')

ERROR:

  ERROR: value "042159162153000033" is out of range for type integer

Below my XML from the REQUEST POST

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"><wfs:Insert><feature:dbtcivici xmlns:feature="http://www.tinyows.org/"><feature:the_geom><gml:MultiPoint xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:900913"><gml:pointMember><gml:Point><gml:pos>1141040.5466402 5682597.1000195</gml:pos></gml:Point></gml:pointMember></gml:MultiPoint></feature:the_geom></feature:dbtcivici></wfs:Insert></wfs:Transaction>

Does anyone know how to solve this problem?


texts: Identifying selected feature using PyQGIS?

texts: If I select a feature and run self.iface.actionIdentify().trigger() it will open the identify window but won't open the selected feature's data, how can I make it open the data of the selected feature?


texts:
I created a webmap with qgis leaflet. It's working well on my pc. I uploaded the folder to my ftp server but on the internet it has no layers, only the osm basemap. I have no idea why.

How can I solve this problem?


texts: Why is the Frequency Tool in ArcGIS so Fast and Efficient?

texts:
Why is the Frequency tool in ArcGIS is so fast and efficient?

Is it possible to make my own Frequency Tool with similar performance (in python)?

I tried to do some based on python dictionary and searchCursor but when I have more than 100,000 records to analyse it takes a long time.


texts: How to calculate polyline azimuth or polyline direction?

texts:
I am looking for a Python code block to calculate azimuth per each polyline feature.
So far I have been able to use the following Python code to calculate bearing:

180  + math.atan2(( !Shape.firstpoint.X! - !Shape.lastpoint.X! ),
(!Shape.firstpoint.Y! - !Shape.lastpoint.Y! ) ) * (180 / math.pi)

However bearing is not sufficient for me as I am looking to ultimately find line directions which can be detected by caclulating an azimuth.


texts: How can I access a specific button in a toolbar using PyQGIS?

texts:
I'm trying to "click" a button in the "Raster Toolbar" in QGIS 2.14 via PyQGIS, but don't know how to do this correctly. I can access the toolbar itself with

rtools = qgis.utils.iface.rasterToolBar()

But after that, I have no clue how to get the correct button. I think I can aceess the single buttons like this:

buttons = [c for c in rtools.children() if type(c) is QToolButton]

but I only know the hover text (e.g. "Local Histogram Stretch") and not the name or ID or whatever can be used here as an identifier. All fields which I thougt might contain something useful have been empty strings.

Does anyone know the proper way of "clicking" a button via code?

Just found the actual action I was searching, so I could use Trigger plugin button using PyQGIS? to get that specific button to be "clicked", but what if I don't know the QAction?

There must be a way of identifying the button through one of its properties.


texts:
currently I'm using PostgreSQL/PostGIS for a small spatial project. Last week the project manager asked me if we could migrate to MySQL with Spatial Extension.

I know MySQL but have no experience with the Spatial Extension. I heard that the extension is not very fast and might contain some bugs.

Now my question - is MySQL Spatial extension stable and fast enough to be considered for an online database with about 20000 spatial operations a day?

The operations are just the standard stuff like:
- export to KML
- point in polygon
- overlapping polygon

The database has one GB of spatial data.


texts: Quickly digitising several polygons based on classified raster output using ArcGIS Spatial Analyst

texts:
I have a raster DTM of a river which I would like to split into 90 sections, each measuring 125m. I have used Euclidean Distance to outline each section by reclassifying the symbology to change colour at each appropriate interval. I would now like to extract each coloured section from the Euclidean Distance symbology to create individual rasters.
Is there a quick way to do this without having to manually cut around each section?
Here is what I am looking at


texts: Looking for a work-around to AGOL's 30-minute cache

texts:
Is there a way to bypass AGOL's 30-minute cache? We are serving up kmls through ArcGIS Online for power outages on an electric network and we need quicker updates than the set 30-minute time frame. Is this possible?

We've read that the KML utility service (http://utility.arcgis.com/sharing/kml) can be tweaked to ignore the cache, but we're not sure how to do this - any ideas?

I have looked for an answer here and posted a question to GeoNet as well, but have yet to find an answer.


texts: Using ArcGIS API interpolation functions (idw, kriging, splines) in my own dataset?

texts:
Does ArcGIS API provide interpolation functions, like if I want to interpolate from my own set of values?

I know that ArcGIS provides functions for interpolation like IDW, Spline and Kriging etc, but can I use those functions to interpolate against my own dataset?


texts: Is it possible to display all icons and labels but not cover other icons in QGIS?

texts:
I've been trying to find a way to automatically distribute icons and labels from different data sources so that they don't cover each other at any scale.
The PDF will have the data showing at this scale and so I can't just reduce the size of the icon or it won't be visible.

Issue 1 - this label should move to the left to not cover the icons.
Issue 2 -these icons should automatically distribute with an arrow rather than be on top of each other.
The settings I have for 1 are basically

I know this can be done with the label move options to an xy position but then when the map zoom level changes this would get affected and we need to print at different scales.
The settings for 2 are basically

I would like this to look like the following which is the layout exported as SVG and then edited in Inkscape

Is there a geometry generator code that would automatically distribute the icons and add an arrow (from the closest vertex) to indicate where it's supposed to be?
Another solution could be if there is a way to just move labels and objects in the print composer without impacting the placement in the map.
Here's an example of what Solution 1 for the icon issue from Corentin Lemaitre looks like


texts: Detect coordinate system UTM 32 North or UTM 33 North

texts:
I am getting database tables (dbase) from customers with geo coordinates either in Gauss-Krueger (zone 2 to 5) and UTM 32 North or UTM33 North, usually from within Germany. Sometimes there is no explicit description which coordinate system is being used (and the customer sometimes doesn't even know it himself) so I would like to detect it.

Gauss-Krueger is easy because the easting value is prefixed with the meridian / zone, so I know that e.g. 4584282;5836214 is Gauss-Krueger Zone 4.

With UTM that's apparently not the case (no preceding zone number 32 / 33 in the data). Is there any way to automatically detect the zone? E.g. how do I know where 381208;5835042 is located? In this case I know for sure that it is zone 33, near Berlin, but in other cases that's not known. Any hints?


texts: Orfeo Toolbox PCA fails - Null Eigen value?

texts:
I am using the Orfeo ToolBox to do a PCA on a .tif with 3060 timesteps of temperature values.

The PCA runs either in the Ubuntu CMD line or in QGIS.

In both, I get the same error saying:

  Null Eigen value

I already found the lines which are throwing the error in the original Python file but I can't figure out what is exactly going wrong.

Can anyone explain this error to me?


texts: PyQGIS error: 'QgsVectorLayer' object has no attribute 'xMinimum'

texts:
I have the following code:

def moveFiles(self):
    mask = self.getBoundingBox() #BoundingBox of polygon layer
    inputDir = self.dlg.input_folder.text() #directory of xyz files - "Users/denni/Desktop/xyz_data"
    outDir = self.dlg.output_folder.text()  #output directory - "C:/Users/denni/Desktop/output/"
    os.chdir(inputDir)
    crs = "EPSG:25832"

    for file in glob.glob("*.xyz"):
        uri = "file:///" + inputDir +"/"+ file + "?type=csv&delimiter=%s&crs=%s&xField=%s&yField=%s" % (" ", crs, "field_1", "field_2")
        name = file.replace('.xyz', '')
        lyr = QgsVectorLayer(uri, name, "delimitedtext")
        outputPath = outDir + file

        if mask.xMinimum() <= lyr.xMinimum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMinimum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions=['GEOMETRY=AS_XYZ'])
        if mask.xMinimum() <= lyr.xMaximum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMinimum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions=['GEOMETRY=AS_XYZ'])
        if mask.xMinimum() <= lyr.xMaximum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMaximum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions=['GEOMETRY=AS_XYZ'])
        if mask.xMinimum() <= lyr.xMinimum() <= mask.xMaximum() and mask.yMinimum() <= lyr.yMaximum() <= mask.yMaximum():
            QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV", layerOptions=['GEOMETRY=AS_XYZ'])

and the following error message:

  'QgsVectorLayer' object has no attribute 'xMinimum'

In the first line, I get the polygon.boundingBoxOfSelected() of a polygon. boundingBoxOfSelected() returns a QgsRectangle according to the documentation . And according to this documentation, I should be able to access the functions xMinimum() and xMaximum().

Thus, I don't understand the error message. What am I doing wrong?


texts:
Im writing a stand alone app in Python using QGIS 2.01 on Ubuntu via Eclipse bindings.

I've written code which successfully adds features to a new in memory layer and displays, but now I want to add attributes. The following code is giving an error : "Undefined Variable: QgsField" error despight importing both "qgis.core *" and specifically QgsField.

import os
import os.path
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.gui import *
from qgis.core import *
from PyQt4.QtGui import QAction, QMainWindow
from PyQt4.QtCore import SIGNAL, Qt
from qgis.core import QgsVectorLayer, QgsField, QgsMapLayerRegistry
from PyQt4.QtCore import QVariant

Graphiclayers =  QgsMapLayerRegistry.instance().mapLayers().values()
        GraphLyr = Graphiclayers[0]

        pr = GraphLyr.dataProvider()
        # add fields
        pr.startediting()
        pr.addAttributes([QqsField('Colour', QVariant.Int),
                          QgsField('QTriad', QVariant.Int),
                          QgsField('QBox', QVariant.Int)])


texts: What SQL query will replace nth character(s) in QGIS DBManager?

texts:
I have a table with a date/time column.  The symbols separating the date need to be changed in order for it to work with the TimeManager plugin.

Eg. 2016:08:13 08:41:20 to 2016/08/13 08:41:20

The column name is 'datetimeor' and the table is 'ba396'


texts: Sum the area of buildings inside a block (QGIS)

texts:
My problem is to get the sum of the area of all buildings inside a block.

There are a layer with the buildings (they contain an block_id to which block they belong and the area) and a layer with all blocks (they have a block_id). I want to sum the areas of all buildings with the same block_id and store this sum in a new column in the block-layer to the block with this block_id.
I add an picture to clarify my question and structure of data. There are two blocks. The left on has the block_id = '078060' and the right one block_id = '078070'. Now all my buildings in the left polygon have an attritbute also called block_id with the value '078060'. I want then to sum up all the areas of the buildings in the block and store this value to my block layer.

I thought about calculation this with the field calculator in QGIS. But I don't know how.

The question "Summing up the area on a polygon by a certain attribute result" is just a first step to solve my problem, because the calculated sum for the categories (block_id) should be added to a new column in the blocks_layer. I was looking for a possibility to automatically calculate the sum and store it to a new column.


texts: how to use gdal_calc to multiply raster by a float variable

texts:
gdal_calc = 'python ' + basedir + '/gdal_calc.py -A' +wind_img+ " -B" + r2+   " --outfile=" + basedir + "/n16_temp_img.tif " + "--calc=" + '"(A * B)" '+ "--type=Float32" + " --overwrite"
os.system(gdal_calc)
n16_temp_img = 'n16_temp_img.tif'

I am trying to multiply a raster by a float variable in python . if i use a number in the operation it works but once I replace the number with a predefined variable it will break .  (wind.img is the raster and r2 is the float variable )


texts: Clipping a TIGER Census Tract by Vector Outline

texts:
So I have recently downloaded the 2011 Tiger Census Tract for the 5 Boroughs of Manhattan which I would like to process for a future visualization.

The Tract Shp file includes water areas, which I would like to clip out using the actual boundaries of each borough.

Or in other words how do I clip the blue tract areas using the red outlines.

Thanks


texts: Converting Lambert Conformal Conic to WGS84 using Proj.4

texts:
I am attempting to transform data from Lambert Conformal Conic lat/long and I am having some issues with the results.

The shapefile I am attemtping to transform can be found here: http://tfr.faa.gov/save_pages/0_8326.shp.zip

Since I am using proj.4 I used GDAL/srs to convert the WKT to a proj.4 init string. In order to do that I had to change the WKT projection from "Lambert_Conformal_Conic" to "Lambert_Conformal_Conic_SP2".

The WKT I used:

PROJCS["North_America_Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_SP2"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-77.03638889],PARAMETER["Standard_Parallel_1",38.69259533705237],PARAMETER["Standard_Parallel_2",39.02625053604863],PARAMETER["Latitude_Of_Origin",40.0],UNIT["Meter",1.0]]

The resulting proj.4 string from the WKT:

"+proj=lcc +lat_1=38.69259533705237 +lat_2=39.02625053604863 +lat_0=40 +lon_0=-77.03638889 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

The proj.4 string I use for WGS84 lat/long:

"+proj=longlat +ellps=WGS84 +no_defs"

When I get a Lambert Conformal Conic point from the shapefile, I transform it to Lat/Long using this line of code:

pj_transform(pj_lcc, pj_latlong, 1, 1, &fX, &fY, NULL );

Where pj_lcc is my Lambert Conformal Conic projection info, pj_latlong is my WGS84 lat/long projection info and fX, fY are the Lambert Conformal Conic point.

An example result from the above transform is:

lcc point: {4842.47, -71287.5} --- long, lat: {-1.34356, 0.686926}

The lat, long should be around 38º51'34"N, 77º02'11"W (center of the shape according to http://tfr.faa.gov/save_pages/detail_0_8326.html )

Any clues as to what I am doing wrong?


texts: Oracle query using 'with' keyword fails in ogr2ogr

texts:
I am trying to return a spatial resultant by throwing a fairly long and complex query at some datasets in an Oracle database, and returning a resultant via ogr2ogr. The goal is to return a resultant in a GeoJSON format after doing a right excluding join - this saves me the need to later run a Python cursor to weed out records with an ID number found in two different datasets.

Here is the Oracle query, tested successfully against an Oracle 12c database (Warning: it's long):

with table1 as (
     select forcov.OPENING_ID as forcov_OP
    from WHSE_FOREST_VEGETATION.RSLT_FOREST_COVER_POLY forcov,
       WHSE_LAND_USE_PLANNING.RMP_LANDSCAPE_UNIT_SVW lu
  where lu.LANDSCAPE_UNIT_NAME = 'Ladybird'
  and UPPER(forcov.SILV_POLYGON_NO) NOT LIKE 'PAS%'and
      UPPER(forcov.SILV_POLYGON_NO)  NOT LIKE '%WTR%' and
     UPPER(forcov.SILV_POLYGON_NO)  NOT LIKE '%WTP%' and
    UPPER(forcov.SILV_POLYGON_NO)  NOT LIKE '%RD%' and
   UPPER(forcov.SILV_POLYGON_NO) NOT LIKE 'RES%' and
  UPPER(forcov.SILV_POLYGON_NO)  NOT LIKE '%ROAD%' and
  UPPER(forcov.SILV_POLYGON_NO)  NOT LIKE '%LANDING%'
  and ((SDO_RELATE(forcov.GEOMETRY,  lu.GEOMETRY, 'mask=inside') = 'TRUE')
  or (SDO_RELATE(forcov.GEOMETRY,  lu.GEOMETRY, 'mask=OVERLAPBDYINTERSECT') = 'TRUE')))  ,

table2 as (
      select vri.OPENING_ID as vri_op, vri.GEOMETRY as GEOM_vri,
      case when LINE_7B_DISTURBANCE_HISTORY like ('L%') and TO_NUMBER(SUBSTR(LINE_7B_DISTURBANCE_HISTORY, 2,2)) <= TO_CHAR(SYSDATE, 'YY')
        then ('2000' + SUBSTR(vri.LINE_7B_DISTURBANCE_HISTORY, 2,2))
      when LINE_7B_DISTURBANCE_HISTORY like ('L%') and TO_NUMBER(SUBSTR(LINE_7B_DISTURBANCE_HISTORY, 2,2)) > TO_CHAR(SYSDATE, 'YY')
        then('1900' + SUBSTR(vri.LINE_7B_DISTURBANCE_HISTORY, 2,2))
      when LINE_7B_DISTURBANCE_HISTORY is null
        then null
      else null end as YEAR_LOGGED

      from WHSE_FOREST_VEGETATION.VEG_COMP_LYR_R1_POLY vri, WHSE_LAND_USE_PLANNING.RMP_LANDSCAPE_UNIT_SVW lu
      where ORG_UNIT_CODE IN ('DSE')
      and lu.LANDSCAPE_UNIT_NAME = 'Ladybird'
      and vri.LINE_7B_DISTURBANCE_HISTORY LIKE 'L%'

      and (
      (SDO_RELATE(vri.GEOMETRY, lu.GEOMETRY, 'mask=OVERLAPBDYINTERSECT') = 'TRUE')
      or (SDO_RELATE(vri.GEOMETRY, lu.GEOMETRY, 'mask=inside') = 'TRUE')  )  )

select table2.vri_OP, table2.GEOM_vri
from table1 right join table2 on table1.forcov_OP = table2.vri_OP
where table1.forcov_OP is null

I think the with clause may be causing a problem, but it's difficult to know; when I test this in the OSGeo4W window using an ogr2ogr command line, it just gives up, without throwing any error message.


texts: Prevent layer from WMSGetFeatureInfo request on openlayers web map

texts:
I use this example code get feature info from my openlayers web map.I add 4 layers to map including raster image as base map.
But i don't need feature info from that raster layer.But when clicking on the map to get information It gives information of all layers.

How I modify above mentioned example code to prevent information from raster layer?

Pls help me,Any idea welcome.Thanks


texts: Setting value, key pair in Dictionary using Google Earth Engine JavaScript API

texts:
A simple script in the Google Earth Engine code editor,
var dict = ee.Dictionary();
dict.set('num', 10);
print(dict);

gives Object (0 properties) as the output
even when i do print(dict.keys()) or print(dict.values()) it throws empty list.
How do I pass a key,value pair to a dictionary?


texts:
I'm wondering if there's an easy way to make a polygon edge in QGIS 2.18 look like an Autocad revcloud (shown in red).

I would prefer an edge that looks a little random over one that is made up of a single arc.


texts: Add raster image (GeoTIFF) over a TMS with Mapbox

texts:
Having some GeoTIFFs, which I would like to display georeferenced over a TMS (in an impossible way to download the GeoTIFF, like another WMS or TMS). Is that possible?
Note: Currently I use Mapbox with react (react-mapbox-gl), I would prefer to solve it that way, but other methods can count as well (e.g. using PostgreSQL or GDAL).


texts: Separating lots of polygons connected by single vertex into individual polygons using ArcMap

texts:
I've been given a shapefile where lots of polygons are connected by a single vertex.
I've included an image which shows when you select 1 polygon it selects lots as they are all connected and overlapping in some places.
How can I separate the polygons in ArcMap?


texts: Excel file . xlsx can't join with the shapefile map layer

texts:
Trying to join a shapefile layer - polygons - (municipalities with a data table with corresponding data sets) an Excel file with the .xlsx extension, BUT as showed in the attached image some error occurs.

Any suggestions what's wrong with this join procedure?


texts: Same Image and Shapefile are not overlapping in ArcGIS for Desktop?

texts:
An image is Digitized in ArcView 3.3 with JPEG extension; But not with Image analysis extension.

When I open that image and shapefile in ArcGIS 10.1 for Desktop then they are not overlaping each other ... why is that?

How can I overlap the image and shapefile?


texts:
I have imported a shapefile into CartoDB (it is the only one I can get my hands on, since it's a small geographic region, non-US) - it includes shapes of all the municipalities in a federal state. It turns out though that the shapefile creates two seperate maps when imported into CartoDB, because of some bureaucratic distinction the state makes in regard to its municipalities.

I would like to have all the shapes of the municipalities in one map and one layer, so I can join sets of data (that normally just include all the municipalities in alphabetic order) easily and also use the info window function for all shapes.

I tried merge tables, but that function seems to be for different types of merges. I tried opening one of the maps and typing out the missing municipalities in new rows in the data window (to trick the column merge) but it appears that the column that has the geospatial info in it is not supported by the merge function.


texts: OpenLayers 3: GeoJSON feature clustering

texts:
I've been trying to use the basic example from the OpenLayers website to create a simple clustering map for my GeoJSON data. However, I'm not getting any luck generating any clusters.

Below is my code:

var baseLayer = new ol.layer.Tile({
    source: source: new ol.source.OSM()
});
var view = new ol.View ({
    center: [19425838, -4409677], // projection is in EPSG: 3857
    zoom: 6,
    maxZoom: 28,
    minZoom: 1
    });
var source = new ol.source.GeoJSON({object: geojson_jobs}); // local .js file passing geoJSON features into a variable
var clustersSource = new ol.source.Cluster({
        distance: 10,
        source: source
    });
var styleCache = {};
var clusters = new ol.layer.Vector({
    source: clustersSource,
    style: function(feature, resolution) {
        var size = feature.get('features').length;
        var style = styleCache[size];
        if (true) {
            style = [
                new ol.style.Style({
                    image: new ol.style.Circle({
                        radius: resolution,
                        stroke: new ol.style.Stroke({
                            color: "#fff"
                        }),
                        fill: new ol.style.Fill({
                            color: '#3399CC'
                        })
                    }),
                    text: new ol.style.Text({
                        text: size.toString(),
                        fill: new ol.style.Fill({
                            color: '#fff'
                        })
                    })
                })];
            styleCache[size] = style;
        }
        return style;
    }
});
var layers = [baseLayer, clusters];
var map = new ol.Map({
        controls: ol.control.defaults().extend([new ol.control.ScaleLine({})]),
        target: 'map',
        renderer: 'canvas',
        view: view,
        layers: layers
    });


texts: geom_rect rectangle disappears when used together with coord_sf in R

texts:
I am plotting a map with ggplot2 and would like to have a rectangle surrounding it. I have tried to create the rectangle using geom_rect function, but whenever I include coord_sf function in the code chunk, the rectangle disappears. Below are the maps I got with the following codes. Classes of zone_map and ecoregion_polygon objects are "sf" and both are read in using st_read function.

ggplot() +
  geom_sf(data = zone_map, color = "black", fill = gray62_transparent) +
  geom_sf(data = ecoregion_polygon, color = "black", fill = burlywood1_transparent) +
  geom_rect(aes(xmin = 20, xmax = 50, ymin = 33, ymax = 47), color = "black", fill = NA)

ggplot() +
  geom_sf(data = zone_map, color = "black", fill = gray62_transparent) +
  geom_sf(data = ecoregion_polygon, color = "black", fill = burlywood1_transparent) +
  geom_rect(aes(xmin = 20, xmax = 50, ymin = 33, ymax = 47), color = "black", fill = NA) +
  coord_sf(crs = "+proj=laea +lat_0=38 +lon_0=35")

As you see in the second map, the geom_rect rectangle disappears upon addition of coord_sf line. How can I solve this issue? I have tried to create another polygon using st_polygon function, but I obtained the same results. I also tried setting both objects' crs as the azimuthal "+proj=laea +lat_0=38 +lon_0=35" projection, which I want to plot my map in, and I got the same results. Below are the st_crs outputs of zone_map and ecoregion_polygon objects.

> st_crs(zone_map); st_crs(ecoregion_polygon)
Coordinate Reference System:
  EPSG: 4326
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"
Coordinate Reference System:
  EPSG: 4326
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"


texts: Getting Light Gray Canvas and other missing basemaps to show in ArcMap

texts:
I don't have Light Gray Canvas, Dark Gray Canvas, and others in my Add Basemap window.
I am using ArcMap 10.7.1

How can I get those added into my add basemap window?
So that it can be like this -


texts:
I have the following function in order to add a sorting value to a layer:

array_find( array_sort( array_agg("field")), "field")

This works fine if I use it through the calculator within the interface (abacus icon) but the same function doesn't work within the field calculator in modeler. The error sais that arrays can't be use it or something like it.

I need to sort a layer or add a rank attribute to a layer from QGIS 3 modeler. I insist: using modeler.

Some screenshots:

By the way, In the model I'm trying to connect an input "Rank" within the function. That not the problem because I already tried using the simple function I wrote above and It worked. The issue here is the "array_" functions that, within the models, just doesn't work.


texts: Publishing 3D scene from ArcGIS Pro 1.2 - Extruded Polygons are not working

texts:
I am trying to publish this map to ArcGIS Online:

These polygons are symbolized by color based on a categorical variable and extruded based on a numerical variable.

When I publish this layer it is treated as 2D:

I have adjusted elevation of the layer to "relative to the ground" and 43,000 m.

Any ideas why these polygons are not extruding properly in the web scene?


texts: Is it possible to reach the Spatial Query plugin functionality in QGIS from a Python script (PyQGIS)? If so, how do I accomplish this? Haven't found any documentation for it, even in the C++ API docs.


texts:
I'm trying to query all the objects within 100m and I want the result to be in GeoJSON not a list of GeoJSON. I'm using PostgreSQL 9.1, PostGIS 2.0. Tried this query:

SELECT ST_AsGeoJSON(geocolumn)
FROM geotable
WHERE ST_DWithin(geocolumn, 'POINT(1000 1000)', 100.0);

but the result is a list of GeoJSON.

Is there a way to convert the result of my query to a single GeoJSON output?


texts: How to get VLRs from a .las file with laspy

texts:
I have a .las point cloud data file that I am trying to get the vlr data from. I have found this link - Update EVLR header info for .LAS file using laspy and have read the docs on Laspy - https://pythonhosted.org/laspy/header.html
It seems like the header vlrs should return a variable length record if it exists.
las_file =laspy.file.File(file.las, mode = "r-")
hdr = las_file.header
hformat = hdr.header_format
for spec in hformat:
    print(spec.name, spec.offs, spec.num)

inVLRs = hdr.vlrs
print(inVLRs)
vlr_count = len(inVLRs)
# OGC Coordinate System WKT
print(hdr.max, hdr.version, vlr_count)

I never get anything back for the vlrs.  I get all other header info without an issue. Sample file here.
When I use WhiteBoxTools - I do get the correct VLR info. - https://pypi.org/project/whitebox/
Am I using the correct method to get vlrs?


texts: Using IDRISI/TerrSet Python API gives ImportError: No module named win32com.client?

texts:
I am trying to use the Python API for TerrSet (formerly Idrisi) to run some tests. According to the downloadable API documentation (https://clarklabs.org/download/idrisi-applications-programming-interface-api-for-terrset/) you only need to have installed TerrSet on your system and run the following:

import win32com.client
IDRISI32 = win32com.client.Dispatch('IDRISI32.IdrisiAPIServer')
...

However, I get an error on the first line of code:

ImportError: No module named win32com.client

Looks like the module cannot be found.

Can anyone help me out on what is happening and how I can access the API in Python?


texts: Downloading specific OpenStreetMap data by tag?

texts:
I'd like to write an iPhone app which uses OSM data to map canals. The data is tagged in OSM as canals, but I'm not entirely sure how to download it. Libraries like route-me don't seem to provide a way to download the data. Ideally I'd like to give the user the option to download certain counties, then scrape certain OSM data tag types.

How would I go about this?


texts: Why won't ggplot fill heat map using geom_tile in R?

texts:
I am trying to emulate some geom_tile plots I saw in a research paper.  The plots were created with simulated data.  Below I created a dataframe of "dummy data" using a subset of the actual simulated data from the research paper (just cut and paste).  And then I do the same thing with my real data.  The real data shows up as very tiny dots in the plot, whereas the dummy data fills the heat map, as expected.
library(ggplot2)

library(gridExtra)

## Create data frame of dummy data and check structure

UTMEasting <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
UTMNorthing <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
sp1 <- c(1, 0, 1, 0, 0, 1, 1, 1, 0)
x1 <- c(-1.69,-1.48,-1.35,-1.44,-1.72,-0.38,-0.64,-0.31,-0.82)
dummy.data <- data.frame(UTMEasting, UTMNorthing, sp1, x1)
str(dummy.data)

## Data structure of dummy data looks fine:

# 'data.frame': 9 obs. of  4 variables:
# $ UTMEasting : num  1 2 3 4 5 6 7 8 9
# $ UTMNorthing: num  1 1 1 2 2 2 3 3 3
# $ sp1        : num  1 0 1 0 0 1 1 1 0
# $ x1         : num  -1.69 -1.48 -1.35 -1.44 -1.72 -0.38 -0.64 -0.31 -0.82

## Plot Dummy Data

p1 <- ggplot(dummy.data, aes(x=UTMEasting, y=UTMNorthing, fill=x1)) + geom_tile() +
    labs(title="x1") + theme(legend.position="none")
p2 <- ggplot(dummy.data, aes(x=UTMEasting, y=UTMNorthing, fill=sp1)) + geom_tile() +
  labs(title="sp1") + theme(legend.position="none")
grid.arrange(grobs=list(p1, p2), ncol=1, main="Dummy Data")

## Remove these data objects and create new dataframe with real data

rm(UTMEasting, UTMNorthing, sp1, x1)
UTMEasting <- c(629385.0, 629215.0, 629489.0, 629886.0, 640709.0, 640686.0, 640607.0, 640515.0, 640369.0)
UTMNorthing <- c(5113783, 5114028, 5113735, 5114063, 5073276, 5073224, 5073390, 5073558, 5073921)
sp1 <- c(1, 1, 0, 1, 1, 1, 1, 1, 1)
x1 <- c(5.586717, 5.517154,5.622111,5.798146,5.814360,5.803089,5.842446,5.878588,6.000669)
real.data <- data.frame(UTMEasting, UTMNorthing, sp1, x1)
str(real.data) # all good

## Data structure of real data looks fine::

# $ UTMEasting : num  629385 629215 629489 629886 640709
# $ UTMNorthing: num  5113783 5114028 5113735 5114063 5073276
# $ sp1        : num  1 1 0 1 1 1 1 1 1
# $ x1         : num  5.59 5.52 5.62 5.8 5.81

## Plot real data

s1 <- ggplot(real.data, aes(x=UTMEasting, y=UTMNorthing, fill=x1)) + geom_tile() +
  labs(title="x1") + theme(legend.position="none")
s2 <- ggplot(real.data, aes(x=UTMEasting, y=UTMNorthing, fill=sp1)) + geom_tile() +
  labs(title="sp1") + theme(legend.position="none")
grid.arrange(grobs=list(s1, s2), ncol=1, main="Real Data")


texts:
In my ModelBuilder toolchain I created a raster and want to create another empty raster wit the same properties. For this purpouse I choose the "Create Raster Dataset". For this tool I need the following properties of my raster: number of bands and the pixel type.

The Number of bands I can extract with "Get Raster Properties" but I can't get the pxiel type with this tool.
Is there a way to get to the pixel type information?


texts:
Can you please help me to correct link rapresentation when the thickness exceed the links.
In the image you can see the value on two consecutive links that generate a bad view.
I think it depends on the thickness linked to the offset, but I am not able to solve it.
With offset=0 I have no problem, but I need to set the offset to put the value on the driving side. Thanks


texts: mapbasic update a column with centroid of object

texts:
I have a table that has a few polygons. I would like to extract the centroid of each polygon and fill up the columns named, Cent_X, and Cent_Y.

The syntax for getting the centroid is CentroidX(obj) and CentroidY(obj). A statement like this doesn't work:

    Update TableA
        Set Cent_X = CentroidX(Obj)

Can anyone please help?


texts:
I have these lines of code:

import tkFileDialog, arcpy, arcpy.da, numpy as np, Tkinter
import arcpy.mapping as mapping

fc = tkFileDialog.askopenfilename(filetypes=[('shapefiles','*.shp')],title='Choose a Shapefile')
print 'fc = {0}'.format(fc)
desc = arcpy.Describe(fc)

#Find the number of rows(alternatives) and print it!
number_of_rows = arcpy.GetCount_management(fc)
R = int(number_of_rows.getOutput(0))
print R
#Find the number of columns(criterias) and print it!
C = len(arcpy.ListFields(fc))
print C

It creates a GUI where i can browse and open a shapefile. This shapefile is an ArcMap shapefile with a layer and an attribute table. With this algorithm i get this error:

Traceback (most recent call last):
  File "C:\Users\Th\Desktop\Python Scripting\Σμαραγδάς Αθανάσιος     Διπλωματική\Python Codes\Python_Codes\Script1.py", line 8, in <module>
    number_of_rows = arcpy.GetCount_management(fc)
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 13613, in GetCount
raise e
ExecuteError: ERROR 000229: Cannot open C:/Users/Th/Desktop/waterways.shp
Failed to execute (GetCount).

How can i solve this? I can access the file with the line below and do the same work with it, but i want to create a GUI so that i can open any file i want to by browsing.

fc = "C:/Users/Th/Desktop/waterways.shp"

I am trying to open the file with the GUI and then use it for any reason i want.But after it prints the directory of the file(as the algorithm demands), i get the error above.


texts: OpenLayers plugin in QGIS, error in projection 3857 Pseudo Mercator

texts:
I'm using QGIS plugin OpenLayers to set the background of Google satellite under my vector layer. However, even if I project my layers as EPSG:3857 (the default and only possible projection for OpenLayers) my layers do not overlap correctly.

What is strange, is that at the beginning it was working, then I saved and closed the project and now I can't overlap them correctly any more. In addition, when it was working correctly, I saved the image and in the resulting picture the layers weren't correctly overlapped...

Any suggestion?


texts: Extracting z values from shapefile using OGR API?

texts:
Is there an example of how we can extract the z value from a shp file given the x,y coordinate value (namely utm) using gdal/ogr c++ lib?

alternatively, would it be easier if we convert the shp file into raster layer format and read off value from there?

Would it be easier?

Is there some examples of how this could be done?


texts:
Probably a very basic question here, I want to change the draw order of a shapefile in QGIS, how can this be done?

As you can see in the picture below, only the "10" layer is shown, the others are below it... However, since this is a depth map of Skagerrak, I would like the other layers to show up on top of it.


texts: Debugging NameError: name 'NbrCicle' is not defined?

texts:
I am trying to run the focal statistics tool on this raster for some reason I am getting this name error, even though the script ran fine a moment earlier.

FocStatCalc = arcpy.sa.FocalStatistics(MajFilCalc2, NbrCicle(3, "CELL"), "MAJORITY")
FocStatCalc1 = r"N:\geog327\python15\Project\Results\FocStatCalc.tif"
FocStatCalc.save(FocStatCalc1)
FocStatCalc2 = arcpy.Raster(r"N:\geog327\python15\Project\Results\FocStatCalc.tif")


texts: I have these variables in Leaflet map .. my difficulty is in visualizing the information using a popup in the Leaflet

texts:
var ProvinciasWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Provincias_Moz',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
}).addTo(map);

var DistritosWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Distritos',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
}).addTo(map);

var PostosAdministraticosWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Moçambique_Distritos_e_Postos',
  // styles:'Provincias',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
}).addTo(map)

;
//I tried with this code, but the popup doesn't come out with information ... and it automatically downloads a file ... I don't know what the error is
map.addEventListener("click", onMapClick);
popup = new L.Popup({ maxWidth: 1000 });
function onMapClick(e) {
  var latlngStr = "(" + e.latlng.lat.toFixed(3) + ", " + e.latlng.lng.toFixed(3) + ")";
  var BBOX = map.getBounds()._southWest.lng + "," +
             map.getBounds()._southWest.lat + "," +
             map.getBounds()._northEast.lng + "," +
             map.getBounds()._northEast.lat;
  var WIDTH = map.getSize().x;
  var HEIGHT = map.getSize().y;
  var X = map.layerPointToContainerPoint(e.layerPoint).x;
  var Y = map.layerPointToContainerPoint(e.layerPoint).y;
  var URL = "http://localhost:8080/geoserver/Mocambique/wms?SERVICE=WMS&VERSION=1.1.0&REQUEST=GetFeatureInfo&LAYERS=" +
    "Mocambique:Provincias_Moz,Mocambique:Distritos,Mocambique:Moçambique_Distritos_e_Postos" +
    "&QUERY_LAYERS=" + "Mocambique:Provincias_Moz,Mocambique:Distritos,Mocambique:Moçambique_Distritos_e_Postos" +
    "&STYLES=&BBOX=" + BBOX +
    "&FEATURE_COUNT=5&HEIGHT=" + HEIGHT + "&WIDTH=" + WIDTH +
    "&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&X=" +
    X + "&Y=" + Y;
  popup.setLatLng(e.latlng);
  popup.setContent(
    "<iframe src='" + URL + "' width='400' height='200'         frameborder='0'></iframe>"
  );
  map.openPopup(popup);
}


texts:
In OpenLayers, when a tile is missing in a TMS layer, a semi-transparent pink rectangle shows in its place. Is there a way to load some default image, like one saying "No data for this region under this zoom" instead?

One possibility is to request the tiles through something like a PHP script or a Java servlet. If a nonexistent tile is requested, this servlet would then send the default tile.

However, I'm currently taking the tiles directly from the filesystem, so I'm looking for a solution to be done only with JS / OpenLayers.


texts: Creating folder and subfolders in empty directory using ArcPy?

texts:
I want to create folder and subfolder in an a specific directory by writing code in python.

It is important to check directory and folder if exist a specific subfolder, avoid duplicate or overwrite.

I create a loop but when directory or folder is empty it does not go to the loop or it gives me error.

I don't know how can create my folders.

Part of my code is like this:

# Set workspace
env.workspace = r"..\..\Test"

# Set local variables
main_folder_path = r"..\..\Test"
main_foldNam1 = "OutPutSpace"

# Execute CreateFolder
folders = os.listdir(main_folder_path)
for folder in folders:
#    if os.path.isdir(os.path.join(main_folder_path, main_foldNam1)) == True:
    if arcpy.Exists(main_foldNam1):
        print main_foldNam1, "folder is Exist", "\n"
    else:
        arcpy.CreateFolder_management(main_folder_path, main_foldNam1)

gdb_name = "CaProduc.gdb"
sub_folder1_path = r"..\..\Test\OutPutSpace"

geoDataBs = os.listdir(sub_folder1_path)
for gdb in geoDataBs:
    if arcpy.Exists(gdb_name):
        print gdb_name, " is Exist", "\n"

    else:
        arcpy.CreateFileGDB_management(sub_folder1_path, gdb_name)
        print gdb_name, "\n"


texts: Any free Tools to display more than 5000 geography objects in SQL Server?

texts:
SQL Server 2008 and up allow you to create a spatial index. If you select records that include a geography or geometry object, a spatial data tab appears that can display up to 5000 points or lines or other shapes. You can select a label column, pan, and zoom. I am trying to debug a road network dataset that has disconnected islands of roads and need to be able to display all the line segments for a given state and label the disconnected ones so I can see where the missing roads are. That means I need to display more than 5000 objects.

Are there any free tools that can do this, including the labeling (color labels would be ok if text labels are not supported)? I have about 250,000 line segments, but anything that can display 50,000 would be okay for my work. A plugin for SSMS would be best, but an external tool would work, too. If there is a C# .NET library, I could use that, too.


texts:
When I browse images from https://earthexplorer.usgs.gov/, I can see that all of the images from Landsat are all tilted/rotated at a certain angle.

Why are they rotated? And what's that for?

Below is an image I got from https://earthexplorer.usgs.gov/.


texts: pyqgis 2.18.16 - Segmentation fault error at 2nd run

texts:
When use pyqgis+python for auto-mapping, I encounter a tricky Segmentation fault (core dumped) error, at the 2nd run of the program. I checked previous similar post but not able to fix the issue.

I have a simple function: test_qgis.py, with content:

def test_QGIS_init():
    app = QgsApplication([], True)
    QgsApplication.setPrefixPath('/usr', True)
    QgsApplication.initQgis()
    QgsMapLayerRegistry.instance().removeAllMapLayers()
    ##do the processing, cropping polygon, produce map, etc.
    ##app.exitQgis()

and then test with a test file "test.py" in the same folder:

import test_qgis
print '-------------------test qgis start -------------'
test_qgis.test_QGIS_init()
print '-------------------test qgis 2nd round -------------'
test_qgis.test_QGIS_init()
exit(0)

The output will be fine for the 1st run, but failed at the 2nd time.

-------------------test qgis start -------------
-------------------test qgis 2nd round -------------
Segmentation fault (core dumped)

Any clue?

It it also interesting that, 1) turn on the exitQgis() at the end of the function does not help; 2)if I run the test through X server remotely, then the test passes without the Segmentation fault error.

I would try to keep this version of pyqgis as some function I have tested only working on this version.

Here is my system information:
OS version

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 16.04.6 LTS
Release:        16.04
Codename:       xenial

python

$ python
Python 2.7.12 (default, Aug 22 2019, 16:36:40)
[GCC 5.4.0 20160609] on linux2

QGIS

>>> from qgis.core import QGis
>>> QGis.QGIS_VERSION
u'2.18.16'
>>> QGis.QGIS_RELEASE_NAME
u'Las Palmas'


texts: Understanding the values from transition layers produced by the R package `gdistance`

texts:
Background and Aim

I would like to create a cumulative raster representing walking time across a surface.

For my purposes, the cost to be accumulated has to be the walking time as influenced by terrain slope, according to the Tobler's hiking function below (which has to be be modified to express time instead of speed):

The corner stone for what I am after is a raster representing terrain elevation (let's call it dtm).

Where I am stuck

I seek some elucidations on the results I got in the very preliminary steps I took to move toward my goal. In particular, I am at loss of understanding the values of what are called transition layers produced by the gdistance package in R.

I am going to provide a reproducible example of what I did (working with a small raster representing terrain elevation), and I will finally point out what I am having difficulties to grasp.

Reproducible example

(1) Load a sample dataset, which will be subsetted to have a smaller raster to work with:

library(gdistance)
dtm <- raster(system.file("external/maungawhau.grd", package="gdistance"))
subR <- crop(dtm, extent(dtm, 30, 36, 50, 56))

(2) Calculate the slope using gdistance; first, the function altDiff defines the difference in elevation between adjacent cells; secondly, the function is applied to the subR raster (elevation) via the transition() command, considering moves in 8 directions; finally, geoCorrection divides the differences by the inter-cells distance, eventually arriving at (if I am not mistaken) slope expressed as rise/run:

altDiff <- function(x){x[2] - x[1]}
hd <- transition(subR, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd,scl=FALSE)
plot(raster(slope), sub="slope (rise over run)", cex.sub=0.8)

(3) The code below calculates the pace (hours per meter) by reworking the Tobler's function and applying the latter to the slope raster; multipling by 1000 and taking the reciprocal helps turning speed (kmh) into pace (hours per meter) (please, disregard the error in the image subtitle, which shows kmh instead of hours per meter).

tobler.hm <- function(x){ 1 / ((6 * exp(-3.5 * abs(x + 0.05))) * 1000) }
adj <- adjacent(x=subR, cells=1:ncell(subR), direction=8)
pace <- slope
pace[adj] <- tobler.kmh(slope[adj])
plot(raster(pace), sub="m per hour", cex.sub=0.8)
plotvalues(raster(pace))

The above would be the cost surface (expressing pace, i.e. meters per hour) that should be accumulated.

Question

My issue at this point is to understand the values stored in the two above rasters (slope and pace) devised using the transion() function. While I understand that, by using the transition() function,

"users  define  a  function f(i,j) to calculate the transition value for each pair of adjacent cells i and j" (quoting from https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf),

I am at loss of understanding how those values (in each respective raster) related to each other? For instance, in the slope raster, how does the 0.291 (upper-left cell) value relate to the 0.102 immediately at its right, and to the 0.249 laying immediately below? The same applies to the values in the pace raster. I need to grasp this (possibly very basic) step in order for me to perform additional calculation to eventually accumulate the pace values from a starting location.


texts: Arc Hydro giving Error HRESLT E_FAIL has been returned from call to Component?

texts:
I have a problem with using Arc Hydro. I tried to follow these steps:

App utilities> Set target locations

Raster Data: C:users\seabell\Deskop\weekdata4\
Vector Data: C:users\seabell\Deskop\weekdata4\untitled.gdb

Fill sinks
Flow direction
Flow accumulation
Stream definition: I gave the number in area
Stream segmentation
Catchment grid deliniation
Catchment polygon processing

I receive an error which is:

  Error HRESLT E_FAIL has been returned from a call to a om Component

Can anyone help me?


texts: Layer to KML GP Service using commas for coordinates instead of points

texts:
I have a model that works fine in ArcMap for converting to KML, but when I publish it as a geoprocessing service, only the first attempt from a web application works. All following attempts return KML files with commas instead of points in their geometries:

First run : -72.5123508930792,46.02718407088821,0 -72.89026165964782,45.9512204096942,0 -73.00197113196214,45.60643088717575,0

Second run : -72,5123508930792,46,02718407088821,0 -72,89026165964782,45,9512204096942,0 -73,00197113196214,45,60643088717575,0

I believe that changing the locale setting (decimal delimiter) of the arcGIS server machine might fix this issue, but changing it permanently is simply not an option. I need to find an alternative.

I have tried exporting the model to python and editing the script to change the locale with python (and then return it to its original locale):

locale.setlocale(locale.LC_ALL,"english")
arcpy.LayerToKML_conversion(sourceLayer, kmzFileName, "1", "false", "DEFAULT", "1024", "96", "CLAMPED_TO_GROUND")
locale.setlocale(locale.LC_ALL,"")

as outlined here: https://community.esri.com/thread/110968

But this doesnt seem to work.

Does anyone know of a way to manually override the way ArcGIS Server/arcpy.LayerToKML_conversion uses the system locale/decimal delimiter?


texts:
I try get number of features through WFS request:

    $.ajax('http://x.x.x.x/geoserver/wfs',{
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typeName: 'workspace:layer',
                outputFormat: 'text/javascript',
                srsname: 'EPSG:4326',
                resultType: 'hits',
                bbox: extent.join(',') + ',EPSG:4326'
                },
            dataType: 'jsonp',
            jsonpCallback:'parseResponse',
            jsonp:'format_options'
    }).then(function(response) {
    var result = parserJSON.readFeatures(response);
    });

But I get Uncaught SyntaxError. Response is not JSONP.

response:

<wfs:FeatureCollection xmlns:xs="http://www.w3.org/2001/XMLSchema"
    xmlns:sf="http://www.openplans.org/spearfish"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:gml="http://www.opengis.net/gml"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:ows="http://www.opengis.net/ows"
    xmlns:tiger="http://www.census.gov" xmlns:topp="http://www.openplans.org/topp"
    xmlns:baikalgeol="http://www.geol.irk.ru/baikalgeol"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" numberOfFeatures="35"
    timeStamp="2016-10-07T02:38:43.015Z"
    xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"/>

As see server return numberOfFeatures="35".

I changed outputFormat on 'GML3' and parseGML3 (openlayers format). But it did not help.

How wrap response in json or read this xml without rest?


texts: Is convention to denote Building Elevation with respect to ground or sea level?

texts: Have a query on shapefile. I have a shapefile containing a number of polygon features representing buildings. I am trying to extract the elevation of the building. I would like to know if the usual convention in denoting the shapefile elevation is given with respect to the sea level or is it given with respect to the ground?


texts: Error in executing clip with polygon from another layer

texts:
I have used clip with polygon from another layer with the Digitizing tools toolbar successfully before but I am running into an error recently.

I have selected one polygon in the layer to clip and one polygon in the clipping layer. Both layers have the same CRS (they are actually polygons copied from the same layer initially just so I could perform the clipping operation).

I do not have openlayers plugin installed (potential issue raised in How to fix a "arguments did not match any overloaded call" error)

  TypeError: QgsCoordinateTransform(): arguments did not match any
  overloaded call:    overload 1: too many arguments    overload 2:
  argument 1 has unexpected type 'int'    overload 3: argument 1 has
  unexpected type 'int'    overload 4: argument 1 has unexpected type
  'int'    overload 5: argument 1 has unexpected type 'int'

  Traceback (most recent call last):   File
  "C:\Users\Kellina\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\DigitizingTools\tools\dtclipper.py",
  line 107, in process
      clipperCRSSrsid, projectCRSSrsid TypeError: QgsCoordinateTransform(): arguments did not match any overloaded call:
  overload 1: too many arguments   overload 2: argument 1 has unexpected
  type 'int'   overload 3: argument 1 has unexpected type 'int'
  overload 4: argument 1 has unexpected type 'int'   overload 5:
  argument 1 has unexpected type 'int'

  Version de Python : 3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12)
  [MSC v.1900 64 bit (AMD64)]  Version de QGIS : 3.0.1-Girona Girona,
  a86bec25eb

  Chemin Python : C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python
  C:/Users/Kellina/AppData/Roaming/QGIS/QGIS3\profiles\default/python
  C:/Users/Kellina/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
  C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins C:\Program Files\QGIS
  3.0\bin\python36.zip C:\PROGRA~1\QGIS3~1.0\apps\Python36\DLLs C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib C:\Program Files\QGIS 3.0\bin
  C:\PROGRA~1\QGIS3~1.0\apps\Python36
  C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\site-packages
  C:/Users/Kellina/AppData/Roaming/QGIS/QGIS3\profiles\default/python
  C:\Users\Kellina\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\DigitizingTools\tools
  J:/MRC Granit/identification des MHH


texts: How to suppress moveend event while switching the base layer?

texts:
I have created a map that contains a two base layers (Google Street Map and Open Street Map). A moveend event listener is registered to listen for any map movement which would perform a function call to load in new data based on the bounds of the map view port. I realise that the function call is performed even when switching the base layer without causing any map movement. Is there a way to switch the base layer silently?

As requested, the code in brief:

var map;
var gmap_r = new OpenLayers.Layer.Google('Google Street Map');
var osm = new OpenLayers.Layer.OSM('Open Street Map');
var options = {
  layers: [gmap_r, osm, markers]
...
};
map = new OpenLayers.Map('map', options);
map.events.register('moveend', map, my_function);
//code above used for triggering a call to my_function upon a map pan or zoom event

The change of base layer unfortunately triggers the above as well.


texts: CssParameter stroke-width min max value Geoserver

texts:
I wonder if there is a minimum value for the stroke width parameter since I've a LineSymbolizer with stroke 0.1 but i want it to be thinner so i change the value to 0.05 but i dont see any difference on the map when the lines are drawn.

PS: sorry my english


texts: GeoServer WMS-layer in OpenLayers only shows labels

texts:
I'm new to GeoServer and encountered the following problem that i don't know how to solve. I have PostGIS-Data in the luxemburghish SRID EPSG:2169 that I transform to EPSG:4326 (seems to be easier to use in OpenLayers) like this:

select f.objectid, f.geocode, f.btyp1_code, f.bewertung1,f.bkjahr,f.shape_ares, st_transform(f.the_geom,4326) as the_geom
from btk_f f

This query is used as a datasource in GeoServer to create a layer, that is served as WMS to be displayed in an OpenLayers Map. The data is styled using an SLD created in QGIS.

In the layer preview in GeoServer the layer displays as expected (showing polygons and labels following the style rules from QGIS)

When I'm loading the layer in OpenLayers (v.5.3.0) only the labels are displayed, not the polygons. Even the labels stop rendering at a certain zoom level.
The resulting map (on OSM as background map) looks like this:

The code I use to insert the WMS-layer follows the latest tiled WMS-example from OpenLayers code.

// only the layer definition part
 var layers = [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    new ol.layer.Tile({
      source: new ol.source.TileWMS({
        url: 'https://ourdomain.com/geoserver/wms',
        params: {'LAYERS': 'store:layer', CRS:'EPSG:4326'},
        serverType: 'geoserver'
      })
    })
  ];


texts: Where to download Tracking Server 10 trial version?

texts: I need to learn tracking server 10. I couldn't find any download link to download trial version of tracking server. Where can i download tracking server 10?


texts: Distances in {Plate Carrée projection - EPSG:4326 - WGS84}

texts:
"4326 is just the EPSG identifier of WGS84" [3]

Indeed http://spatialreference.org/ref/epsg/4326/html/ says it is a GEOGraphic Coordinate System:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

Citing [1], [...]

  one of the most common ways the round world is displayed on a map is
  using the simplest projection we have:

x = longitude
y = latitude

  The name of this projection is “Plate Carrée”, and is widely used
  because it is so simple.  However we often seem to forget that we are
  talking about a projection. Therefore the spatial reference for this
  projection is very often (mis)referenced as a spherical coordinate
  system like for EPSG:4326.

In [1] then a correct WKT for the Plate Carée projection is reported (a projection is overlaid over WGS84/EPSG:4326):

PROJCS["World_Plate_Carree_Degrees",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],
    PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
    PROJECTION["Plate_Carree"],PARAMETER["False_Easting",0],
    PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],
    UNIT["Degree",111319.49079327357264771338267056]]

which is ~EPSG:54001, except that units are not meters but degrees.

I am quite a bit confused for example if, at least in practice, EPSG:4326 is equivalent to Plate Carrée equirectangular projection.

E.g. in a computing environment, projected data is much more convenient since distances are Euclidean and not great-circles: would the distance between two points in EPSG:4326 be different then with respect to the same points in a Plate Carrée?

This is maybe specific to the single application, e.g. [2] Open Layers uses the term 'EPSG:4326' to mean the Plate Carrée projection. R seems to assume a spherical WGS84 instead (-> r-sig-geo).

References:

[1] http://www.sharpgis.net/post/2009/02/06/Why-EPSG4326-is-usually-the-wrong-e2809cprojectione2809d.aspx
[2] Projecting EPSG:4326 data in 2D map?
[3] Difference between WGS84 and EPSG:4326?


texts: Creating clusters based on nearest distance attribute in ArcMap using ModelBuilder

texts:
I am attempting to create a tool via modelbuilder that creates clustering for a polygon feature class based on the distance field attribute (between one polygon to another) that I have created using the "near table" tool. There are some requirements to the clustering process:

Firstly, I would like the cluster to be grouped in an "ascending" method, where the "first" polygons to be grouped together are those with the shortest distance. (i.e. Polygon A and Polygon B with distance of 200m (shortest distance) are grouped together first)

Secondly, I would also like the grouping to be such that once a set of polygons are grouped together in a cluster, they are excluded from the clustering process, which continues on until all polygons are grouped in their own respective clusters based on the distance field attribute.

I would like to know if there is are methods within Arcmap that ultimately allows me to create the clusters through the use of a tool created using modelbuilder, as the existing geoprocessing tools for creating clusters do not have the options that allows me to fulfill the 2 requirements above.
My arcgis desktop is advanced license.


texts:
I having difficulty with the wedge buffering tool in qgis (Wedge Buffer Tool). I applied the tool to a point shapefile (with about 15 points) that I created in arcgis pro.

Problem #1 is that the wedge buffers are all on top of each other - the vertices are based on only one point, whereas I assumed they would be based on each of the 15 points.

Problem #2 is that they are not mapped in the correct location. The original 15 points are located in Massachusetts and the wedge buffers are located in the midwest. They seem to be in the same projection, but I am very inexperienced in qgis. I use arcmap everyday, but decided to try qgis for the wedge buffering (which doesn't seem to be an option in arcgis unfortunately).

I can provide shapefiles if helpful.


texts: Returning a promise from a function in ArcGIS API JavaScript version 4.x

texts:
How would I return a promise from a queryTask from a function in ArcGIS API version 4.x?

function getPromise() {
  require(["esri/tasks/support/Query", "esri/tasks/QueryTask"], function (Query, QueryTask) {
    var query = new Query...
    var queryTask = new QueryTask....
    return queryTask.execute(query)
  });
}

getPromise().then(function(response) {
    .......
    // getPromise() returns undefined here
})

the above would cause an error :

Uncaught TypeError: Cannot read property 'then' of undefined.

I noticed that if I place the require.... outside of the getPromise function, then the promise returns fine. However, when the require... is placed inside the function, it causes the error.


texts:
I am using an application called Gaia GPS to analyze regions affected by large enterprises. So, when I travel to this specific places I mark waypoints and take pictures before, during and at the end of the operations. With this I can analyse lots of factors that could affect the socioenvironmental areas.
Gaia Gps gives me a KMZ file with everything inside, waypoints, routes, photos. At this point ok. I could open this files perfectly on Google Erath Pro. When I export this KMZ using Google Earth it gives me all photos with a KML file. The problem is that the pictures doesn't have EXIF geotag info, neither time where the photos where taken. At the end all I have is a bunch of photos without individual georeferenced information.
<Placemark>
            <name>Minha Foto - 12/02/2021 09:57:23</name>
            <description>
<![CDATA[Location: -24,93659, -51,38044

Elevation: 838,0 m

Time Created: 12/02/2021 09:57:23 BRT

]]>
            </description>

    <Point>
        <coordinates>-51.38044207361656,-24.93658958005247,838</coordinates>
    </Point>
</Placemark>

If you pay attention, the excerpt of KML above show that all the information I need is put inside a custom field CDATA, All I need is to record this on the pictures.
Do you know a way to automatically do this (I have lots of pictures)?
I'm trying to solve this without success for a time now.


texts: FeatureClassToNumPyArray/other ArcMap functions in QGIS?

texts:
ArcGIS user here. Since we're all working from home now, ArcMap only works 2 VPN's deep and is basically garbage to us now at our office, so I'm going all in with QGIS.

I would like to know, what is the equivalent of ArcPy's FeatureClassToNumPyArray which does what the name suggests (grabs attribute table fields, length, area, x , y, etc and loads it to a structured numpy array) and is there any good resource to translate to QGIS and it's Python dialect for ArcMap/ArcPy users?

QGIS version 3.12, Python version 3.7.6


texts: Downloading a huge GeoJSON file for features

texts:
I working on an ASP MVC web app that calls a web service to get the polygons (as GeoJSON files) to display on the client side in OpenLayers. These polygons contain in total 47191 points !! Which makes roughly 5Mb of data...

Are there any strategies to transfer big amounts of GIS data through a web service without the size overhead ?

I've tried GZipping the requests but it didn't help much.

Is it possible to transform the GeoJSON files to some binary format to be downloaded and read by OpenLayers ?


texts: Raster Calculator making incorrect calculations in ArcGIS Desktop?

texts:
I am trying to use raster calculator (ArcMap 10.2.2.3552) to apply a poisson generalized linear model to estimate species abundances based on elevation. The outputs don't make any sense to me, and I'm hoping someone can clarify why.

The best fit model (generated elsewhere) is:

Exp( 3.394 + ( 4.717 * "elev" ) + ( -2.602 * ("elev" ^ 2) ) )

Input raster below:

The output, however, doesn't reflect the math I input, and I am wondering if I have screwed up the syntax, or am missing something basic about raster calculator (which works fine for other models I have used on this dataset).  The "elevmllw" raster ranges from -84 to +5 in values (meters) and the output ranges from 0 (which is correct) to 3.49972e+011 (I'm not sure where that number comes from, solving my model equation for that solution gives an imaginary number).

Output:

When trying to figure out what was going on, I poked around noticed all the values made no sense.  Using identify, it's clear that the calculated values for given cells (from "poisson") do not align with what I am trying to calculate for the elevation value from that cell, so it seems like I must be missing something about how to specify this equation.

The calculation runs without any error messages and the environments are shown below (sorry I'm not sure which are most helpful in troubleshooting, so I showed them all)


texts: Understanding MODIS pixel size, CRS projection SR-ORG:6974 in Google Earth Engine

texts:
I am calculating the area of an ice cover by summing the number of pixels having a value above a certain threshold then multiplying by the pixel area. This is done for a collection of images over a number of years.

The resolution of the MODIS dataset I'm using is 500 m. When I use the function ee.Reducer.sum() I set the scale to this value. When my area is calculated however, the pixel size appears to be 463.3 m with areas reflecting this pixel size (0.216 km2 instead of 0.25 km2).

I am not a remote sensing or GIS scientist by trade and am struggling to understand the relationship between projection (when a pixel is inspected it shows, SR-ORG:6974) and scale. This MODIS data set uses a sinusoidal projection which I understand to be area-preserving.  That being said when I look at pixels nearer to the equator, they also appear to have a size of 463.3.

My bottomline is that I need to present the most accurate measure of area and make sure I have not caused an error.

Are the areas I've calculated using this approach correct or should they be based on a 500 m pixel?

My code is quite long and may confuse the main question, but I've included the main parts below to show my general method. It's not run-able but I can provide the rest if need be.

//Create a number of masks and other functions to clean up images and highlight ice, there is no //use of scale up to this point

var maskedModis_ice_AllQAcloud= sample
  .map(maskLand)
  .map(addice)
  .map(getclouddata)
  .map(getcirrusdata)
  .map(getQAdata)
  .map(maskAllQAcloud);

//Function to create threshold and multiply by pixel area to determine ice area

var icethreshold= 3000;
var icefunction = function(image){
  var icejunk= image.select(['Ice_index']);
  var wateretc= icejunk.lt(icethreshold);
  var solidice = icejunk.gte(icethreshold);
  //image = image.updateMask(solidice);
  var area2= ee.Image.pixelArea();
  var iceArea1 = solidice.multiply(area2).rename('iceArea1');
  var iceArea= iceArea1.divide(1000*1000).rename('iceArea');
  return image.addBands(iceArea);

//highlight only pixels marked as ice through thresholding and combine with all filtering steps

var filteredAlliceQAcloud=maskedModis_ice_AllQAcloud.map(icefunction);
print (filteredAlliceQAcloud);

//Apply to specific regions and plot

var chart = ui.Chart.image.series(filteredAlliceQAcloud.select('iceArea'), regions, ee.Reducer.sum(),500)
    .setChartType('LineChart')
    .setOptions({
      //title: 'Beisfjord',
      hAxis: {title: 'Date'},
      vAxis: {title: 'Ice Area, QA and cloud filtered (m^2)'},
      fontSize: 20,
      legend: {position: 'right'}
    });


texts: How to draw a longest and shortest straight lines through centroid in polygon?

texts:
I have a task I thought some of you might have idea how to solve it.

I have shape and a centroid in it. I need to find longest and the shortest straigt line from the borders of the polygon.

I use QGIS, but have access to ArcMap as well.


texts: ASCII Grid file doesnt have data matching NCOLS

texts:
I have an ASCII Grid file that states the following in its header:

ncols 2205
nrows 653
.
.
.
nodata_value -9999

But when going through the rows, the max number of cols is 12 and the number of rows is around 83,000. It seems like the rows are wrapping around but I am not sure why.

ArcMap ingests correctly, but I am trying to understand what is going on in the file.

Has anybody seen this before?


texts: OpenLayers Drag&Drop a pin to select position?

texts:
Is there an OpenLayers example that allows the user to set a position with drag&drop.
I would like to have a pin and set the position after the user drags & drops the pin to a new position.


texts: Building query tool for spatialite data in QGIS?

texts:
I am trying to build some tool for a collegue of mine to query data from a spatialite database while clicking on vector features.

My first guess is, that I could do this with an "action" I define for the layer. The example actions provided for layers unfortunately do not give me any hint how to perform, and more detailed documentation on actions is quite absent.

I need a query, cause each feature in layer has many corresponding entries in the database table (one-to-many).

While performing the action, I would like to pop up simple table which lists the corresponding data sharing the same id as the feature.

In the example action, i can query an attribute of the layer with this:

echo "[% "id" %]"

which then looks like this:

Ok, that's a start.

But what i need is a table which lists all items in a spatialite table with the same id-attribute as the feature.

My data structure:

Layer attribute:
id

Data table:
id | name | year

The layer itself is stored in the spatialite database too, if thats of any help.

Can anyone give me a hint on how to proceed?


texts: Auto vectorization of raster in QGIS 3.14

texts:
Is it possible to do an auto digitize of a raster (downloaded and WMS)?
I would like to create vector polygons out of it.


texts: Conversion of Landsat ToA reflectance to BoA reflectance with an equation

texts:
I am using temporal datasets from Landsat 5, 7 and 8.
Since I have many datasets I would like to include the conversion from ToA(top of atmosphere) to BoA(bottom of atmosphere) reflectance into a script using Arcpy.

Is there an equation to convert the ToA  reflectance to BoA  reflectance?
I have checked the Landsat data users handbook but I could not find any formula:


texts: How to extract raster from .gdb instead of empty polygons?

texts:
I have met with a problem when importing raster data of light intensity. Concretely, the following file contains raster data of Chinese nightlight intensity. I write the following command to import it:
file_temp = st_read(dsn = 'path/tif.gdb', layer = 'noflare_20131')

light = raster(file_temp)

plot(light)

link for replication data: https://www.dropbox.com/s/o8rn7cnb4a5rc6j/tif.gdb.zip?dl=0
However, it seems that I have only extracted an empty polygon into R without values of light intensity. In contrast, when importing the data into ArcGIS, it could be displayed successfully.
I guess that I should add some options to the command st_read() or raster() while have not figured out a specific solution.


texts: How to get the geometry of a polyline to a specific precision?

texts: In python, I use point = geom.asPolyline() but it looks like I am not getting any precision just whole numbers?  How can I get the geometry and a certain precision?


texts:
I want to know what kind of format (syntax) a geospatial rdf should have so I can load it into virtuoso (open source) and do a GeoSPARQL query over it. I found a very small and simple geospatial rdf and wanted to load it into VOS, but it gave me an error as:

  Unable to load RDF graph  from  with
  Content-Type 'application/rdf+xml': 42000: RDFGE: RDF box with a
  geometry RDF type and a non-geometry content

Edit: I uploaded RDF without any instances and VOS loaded it successfuly. But when I add the instance into it, the error pops. Here the whole RDF i use is:

<?xml version="1.0"?>
<rdf:RDF
    xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
    xmlns:rdfs="http://www.w3.org/2000/01/rdf-schema#"
    xmlns:owl="http://www.w3.org/2002/07/owl#"
    xmlns:sf="http://www.opengis.net/ont/sf#"
    xmlns:geo="http://www.opengis.net/ont/geosparql#"
    xmlns:my="http://example.org/ApplicationSchema#">
<rdfs:Class rdf:about="http://example.org/ApplicationSchema#PlaceOfInterest">
    <rdfs:subClassOf rdf:resource= "http://www.opengis.net/ont/geosparql#Feature"/>
</rdfs:Class>
<rdf:Property rdf:about="http://example.org/ApplicationSchema#hasExactGeometry">
    <rdfs:subPropertyOf rdf:resource="http://www.opengis.net/ont/geosparql#hasGeometry"/>
    <rdfs:subPropertyOf rdf:resource="http://www.opengis.net/ont/geosparql#hasDefaultGeometry"/>
</rdf:Property>
<rdf:Property rdf:about="http://example.org/ApplicationSchema#hasPointGeometry">
    <rdfs:subPropertyOf rdf:resource="http://www.opengis.net/ont/geosparql#hasGeometry"/>
</rdf:Property>
<!-- Instance -->
<my:PlaceOfInterest rdf:about="http://example.org/ApplicationSchema#F">
    <my:hasExactGeometry rdf:resource="http://example.org/ApplicationSchema#FExactGeom"/>
</my:PlaceOfInterest>
<sf:Point rdf:about="http://example.org/ApplicationSchema#FExactGeom">
    <geo:asWKT rdf:datatype= "http://www.opengis.net/ont/geosparql#wktLiteral">
        <![CDATA[ <http://www.opengis.net/def/crs/OGC/1.3/CRS84> Point(-83.4 34.4) ]]>
    </geo:asWKT>
</sf:Point>

</rdf:RDF>

VOS Version is 07.20.3217. What's wrong?