texts: Calculation appropriate angular step depth in depthmapX

texts:
I am a new depthmapX user in which trying to calculate segment step depth from specific segment to store cumulative angles for each segment. Is there anyone who can provide advice to implement it?

depthmapX
http://www.spacesyntax.net/software/

Network conditions are shown below:

Network is composed of center lines produced from polygon of alleyway in a city bock
Alleyway has 5 entrances from sidewalk at the fringe of the block and all alleyways are connected
Center lines are produced from Voronoi polygons, so, there are three characteristics in the network: 1) each segment has very short distance (spreading from less than 1 meter to  3 meters) to replicate  shapes of the alleyway, 2) there are a number of segments between actual intersections due to previous characteristics, and 3) due to the nature of Voronoi polygons, each segment is angled slightly though it is actually straight alleyway.

Followings are the manipulations that I did:

Convert dxf file to axial map and fix disconnection of some lines which can be found by angular step depth analysis specifying a link
Convert axial map into segment map and conduct angular step depth analysis
*I confirmed both analyses shows the same result (axial map also can conduct angular step depth analysis)

The results show following issues:

As there are a number of segments between actual intersections with small angle, cumulative angles from selected segment increase gradually
Due to above result, segment angle after turning right or left does not show significant difference before the turning
Expected results are that segments between intersection have the same or similar value and certain amount of cumulative angels differ from the ones before turning

To solve the issues above, following countermeasures have taken:
- Simplified network, for instance, one segment between two intersections has built to confirm calculation result
- However, at T-shape intersection crossing at right angle for instance, cumulative angles of segments after taking turn does not show 1 point increase compared with the one before the turning
- I thought Tulip analysis can be applied to neglect such small angles as it can approximate angles, however it does not work for the analysis for selected segment

The first issues may be caused by the network specifications and the latter may be more general issue.
I wish if someone has solved similar issue before or provide solution or hints to solve it.


texts: Python: Cannot construct empty GeoDataFrame with columns

texts:
My script
import geopandas as gp
test = gp.GeoDataFrame(columns=['geometry', 'costs'])

the error:
File "/usr/local/lib/python3.7/dist-packages/geopandas/geodataframe.py", line 89, in__init__ super(GeoDataFrame, self).__init__(*args, **kwargs)
file "/usr/lib/python3/dist-packages/pandas/core/frame.py", line 348, in __init__
mgr = self._init_dict(data, index, columns, dtype=dtype)
File "/usr/lib/python3/dist-packages/pandas/core/frame.py", line 451, in _init_dict                                                                                nan_dtype)
File "/usr/lib/python3/dist-packages/pandas/core/dtypes/cast.py", line 1196, in construct_1d_arraylike_from_scalar
dtype = dtype.dtype
AttributeError: type object 'object' has no attribute 'dtype'

What could be the reason for that? It always worked for me like that. Now I had a software update and it suddenly does not work anymore.


texts: How to correct the deviation of reprojected MODIS L1B by using the gdalwarp?

texts:
I am working on reprojecting the MODIS L1B data which downloaded from the NASA EOSDIS Reverb website. The data is MOD021KM. I applied the gdalwarp command to convert the projection of the image to the WGS84.

gdalwarp -of GTIFF -r near -t_srs EPSG:4326 -overwrite -tr 0.018728541 -0.0089737573 HDF4_EOS:EOS_SWATH:"MOD021KM.A2013145.1950.006.2014225011320.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB MOD021KM2013145_2.tif

However, the reporjected image by gdalwarp cannot be matched to the Lake mask in figure 1. In figure 2, the image which was projected by ENVI conversion toolkit is identical to the lake mask.

The right-bottom of the two figures shows the image extent via the coordinates of the four corners. It indicates the longitude is error in the figure 1. How do I correct the deviation of the image using gdalwarp.


texts: How to get azimuth and altitude from a 3D surface in PostGIS?

texts:
I would like to know the azimuth and tilt of a polygon / surface in a PostGIS database. The azimuth should be given e.g. from 0° to 360° degrees. The altitude should be between 0° to 90°. An example result table:

sid | azimuth | tilt | geom
---------------------------
 0  |     0   |   0  | ...
 1  |    90   |  45  | ...

So sid 0 would represent a north facing wall and sid 1 a East facing roof with a 45° degree slope.

Update:
I guess the typical approach would look as follows:

Find the surface's centroid
Project a point North of it
Measure the azimuth and angle between the centroid and the projected point

After giving it some more thought and taking into account that all of the involved functions (ST_Centroid, ST_Project, and ST_Azimuth) are not 3D capable, I feel like this may currently not be possible in PostGIS.


texts:
I am working with a large (30GB) 3 band geotiff. I am using the following code to attempt to import the bands into python

import gdal
import numpy as np

raster_path = "C:/Path/To/Image.tif"
raster_dataset = gdal.OpenEx(raster_path, gdal.GA_ReadOnly)
geo_transform = raster_dataset.GetGeoTransform()
proj = raster_dataset.GetProjectionRef()

bands_data = []
for b in range(1, raster_dataset.RasterCount+1):
    band = raster_dataset.GetRasterBand(b)
    bands_data.append(band.ReadAsArray())

However, when I try to iterate over the bands, all my memory gets used and it won't load for hours.

Any suggestions regarding other ways to handle large spatial data in python? I have used ArcGIS and R in the past, but am relatively new to python. I would prefer to work from command line, if possible.


texts: Create an outline based on density of raster

texts:
I have a classified aerial image of small vegetation that is separated enough by ground that it is not overlapping each other. After the classification process (dzetsaka), I removed all the other classes to only show the wanted vegetation. I am wondering if there is a way to draw a "circle" that best outlines the vegetation. This is currently a raster file, but I can convert it to a vector if need be. The random pieces around the main veg should be ground points, but were classified incorrectly. Likewise, the upper right additional piece should be shadow points.
I am unaware if there are any good models out there for this? I was thinking perhaps it could be done in excel? Or if there is a QGIS plugin? Unfortunately, I think I am going to have non-veg points classified as veg, so I am looking for something that will identify and draw the "circle" based on the dense part.
Any suggestions on where I might find some more information on something like this?

------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------


texts: First time using ArcGIS JavaScript API - where do I find the JavaScript API to include (.js) into my solution?

texts:
I'm writing a web service that will interface with ArcGIS.  One of the requirements is to put:

esri.config.defaults.io.proxyUrl = "proxy.ashx";

in the JavaScript portion of the page.  Of course, it says "esri is undefined" if I don't include the library that defines "esri".  Question: Where do I get this library at?  I've searched all over the internet and couldn't find anything!  Thanks!


texts: Why is the background of my Geoserver WMS+Leaflet example grey instead of transparent?

texts:
I am trying to add a WMS layer to my leaflet map. When I add the layer it shows up but the the outside of the layer is grey instead of a world map.

This is my code:

<!DOCTYPE html -->

<title> leaflet webviewer with Geoserver</title>

<link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.1/dist/leaflet.css" integrity="sha512-Rksm5RenBEKSKFjgI3a41vrjkw4EVPlJ3+OiI65vTjIdo9brlAacEuKOiQ5OFh7cOI1bkDwLqdLw3Zg0cRJAAQ==" crossorigin=""/>

This is the result of the code:

I do not understand what I am doing wrong and why the outside of the layer is not transparent.

Any help/ tips are welcome.


texts: How to use all my data from Google Earth on QGIS?

texts: So my problem is this : I been working in Google Earth for 4 years, and always lacked something, last week I found out about QGIS and started experimenting with it. It have much more tools than GE and I wish to transfer all my old progress from Google Earth to my QGIS Version 3.2 Bonn. Is there any way that I can import all my polygons(fields), lines(rivers,roads) into QGIS so I can continue my work with QGIS?


texts: How to calculate spatial accuracy (%) of the same area of two polygon?

texts:
I run land use modelling for study area in 2015 (in ASCII format and convert to polygon shapefile). I also need to compare the result from the model and actual land use in 2015 (percent accuracy of overall land use, land use type A for example).

Would it be possible to do that in ArcGIS 10.1?


texts: Mosaicking many tif files in ArcMap for ENVI

texts:
my question is how to go about mosaicking ~120GB, 1200+ tif files into one raster within ArcMap. I am currently splitting up the region in smaller sets and then plan to mosaic all these together to form the one final file at the end. I'm using the Mosaic to New Raster tool. The reason I need a single file is so I can use ENVI classification functions on it.

Does anyone know if ENVI can use an ESRI Mosaic Dataset? Is there a better way to go about this (e.g. file type, etc)? I'm using unsigned 8bit data, 4 band imagery. It's taking a long time and I don't know if ArcMap will even be able to mosaic the total, final file at the end of this. I feel like this is one of those "I don't know what I don't know" problems. What am I missing?


texts: Getting OSM data for a particular jurisdiction

texts:
I'm trying to make a street map of New York City (and only New York City). I've downloaded OSM data for the area (from here), but it contains lots of data from the suburbs.

Currently, I've loaded this data into PostGIS (with osm2pgsql), but there's nothing obvious which would let me distinguish roads in the city from those in the suburbs. What's the best way to get streets only within the city?


texts: Qgis python ComboBox: Combination of fixed and free Items

texts:
I want to built a ComboBox to select a given plotting scale out of a list or ask for a different scale (e.g. '1:625') as Textinput in the same ComboBox. The combobox is changed to 'editable' and I'm able to fill in the digits. But they are not read.
The list will work fine -  but how must the code looks like to insert the free choosen scale (last item)?

    self.dlg.comboBoxScale.clear()
    self.dlg.comboBoxScale.addItem('1:500', 500)
    self.dlg.comboBoxScale.addItem('1:1000', 1000)
    self.dlg.comboBoxScale.addItem('1:1500', 1500)
    self.dlg.comboBoxScale.addItem('1:2000', 2000)
    ...
    # Item with free digits
    self.dlg.comboBoxScale.addItem('1:')    #????

Any Idea?


texts: How to find points within a polygon in PostGIS

texts:
I have a series of lat and long coordinates I've loaded into a spacial table (Nov01) in PostGIS, and I've created a geometry column in the table. I also have a shapefile (greatermanchester) I've loaded in with the POSTGIS 2.0 shapefile loader. I've double checked that they both have the same SRID of 4326, and then I ran the code

select a.*, b.*
from "Nov01" a
join greatermanchester b
on ST_Intersects(a.the_geom, b.geom)
limit 10

However it doesn't return any results.
My points table looks like:

What might I have gotten wrong?

EDIT:
There definitely should be points within the polygon. I loaded the original csv and shapefile into QGIS and got this:

I also ran explain analyze queries with these results:

EXPLAIN ANALYZE VERBOSE select count(*)
from "Nov01" a
join greatermanchester b
on ST_Intersects(a.the_geom, b.geom)

result:


texts:
I tried to make a python based CGI script that will do some geoprocessing tasks on the server and then return the result to the user after the processing its finished. To do the geoprocessing tasks, I'm using QGIS python library (PyQGIS).
When I developed the code, I manage to create a standalone python script that running as expected,  But when I used the same script published as CGI, its always failed.
For example, below is a simple code that will load raster layer as QgsRasterLayer. In the last part, it will check whether the layer object is valid or not. When I run this as standalone script, it works fine and it will show that the layer is valid (isValid() method return true), but when I run this script as CGI the layer is always not valid (isValidi() return false).
Is there anything wrong with my code? I'm using Ubuntu 20.04.2 LTS, QGIS 3.18, Apache2.
#!/usr/bin/python3
import os
import sys
os.system('export  QT_QPA_PLATFORM=offscreen')
os.environ["QT_QPA_PLATFORM"] = "offscreen"

# Import and initialize Processing framework
from qgis.analysis import QgsNativeAlgorithms
from qgis.core import *
workdir='/path/of/directory/'
ipRasterSin=os.path.join(workdir,'dataset.tif')

# Initializing
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.setThemeName("default")
qgis = QgsApplication([], False)
qgis.initQgis()
qgis.processingRegistry().addProvider(QgsNativeAlgorithms())

sys.path.append('/usr/share/qgis/python')
sys.path.append('/usr/share/qgis/python/plugins')
sys.path.append('/usr/share/qgis/python/plugins/processing')
import processing
from processing.core.Processing import Processing
Processing.initialize()

rasterLayer=QgsRasterLayer(ipRasterSin,"RasterData")
# To check whether the layer is valid or not.
sys.stderr.write('\nRaster = '+str(rasterLayer.isValid()))
qgis.exitQgis()


texts: How to pull a shapefile into R from PostGIS database (in Windows)?

texts:
I'm currently working on pulling data from a PostGIS database for the first time to map in R and need some help. The database contains a shapefile of countries/regions of the world. I am able to pull the attributes into R as a table through the dbReadTable function from the package RPostgreSQL. This is great for manipulating the data (adding, merging, subsetting, etc.), but ultimately, I still need to be able to actually map this bad boy, so I need the data to be spatial in some way (side note: I plan on mapping with ggplot2, so my data needs to be suitable for that). Here's the tiny bit of code I've successfully used to pull in the data as a table/data frame:

#load in required libraries
library(RPostgreSQL)

#connect to database
drv <- dbDriver("PostgreSQL")
DB <- dbConnect(drv, dbname = "dbname", user = "user",
                      host = "host",
                      password = "password", port = "5432")
#pull in shapefile data
world.shp <- dbReadTable(DB, "ne_10m_admin_1_states_provinces")

Again, I have no problems connecting to the database. That said, I need to pull the data in as a shapefile rather than as a data frame. This question seems to be what I'm looking for, but the accepted/top solution confuses me; specifically, I don't know how to tackle this bit of code to pull in the shapefile from the database:

require(rgdal)
dsn="PG:dbname='gis'"

Similarly, I've found various solutions for Ubuntu, but only this one for Windows, which spits out errors for me.


texts: Site viewshed with DSM: removing "false" visibility readings from trees and buildings

texts:
I've been working on modeling a proposed building and running a visibility analysis to determine from what areas the building could potentially be seen using ArcGIS Desktop with 3D Analyst. Using a DSM raster derived from LiDAR (includes all points - ground, trees, buildings, etc.) as the base surface information so that above-ground sight obstructions are accounted for.

My process to find the areas that the new building could be seen from has been to set the proposed building roof corners and peaks as 3D points and use them as "observers" - if the roof corner can "see" an area, then the building can be seen from that area. I used the Visibility tool in the 3D Analyst toolbox and Model Builder to run through the observer points at the building roof (16 all together) to generate visibility rasters and then mosaic them together.

The issue is that using the DSM is a double-edged sword - using it to account for above-ground obstructions has also created speckled areas of "visibility" on treetops and building tops where the new building could be seen from if the observer's eyes were at that height. But obviously, an observer would not be standing on a rooftop or treetop, so these areas are erroneous. See image below (hillshade with red areas being visible from the proposed building site)

How can I remove these areas, besides manual cleanup? The LiDAR is not classified for vegetation and buildings, only ground and water, so using those LiDAR points as a mask to delete the erroneous areas would require some classification work.


texts: QGIS API : Split a first polyline layer by a second points layer

texts:
Is QGIS or the API QGIS able to split a first polyline shape layer by a second points shape layer as show in the figure bellow.

Is this possible?

The goal is to calculate segments length between points.
The points are positioned exactly on the polyline.
In this example, there would be 3 segments.


texts: Replacing values of raster in R to produce landuse-like raster map?

texts:
I have a raster Z with elevation values (0, 0.12, 2.12), which I created myself in R. I want to produce a land-use raster from these elevation values, i.e. if Z<0.12, I assign land-use class 1 to it, if Z== 0.12, I assign land-use class 2 and if z>0.12, land-use class 3 and so on. I realized that using the method

Z[Z<0.12]<-1,
Z[Z==0.12]<-2,
Z[Z>0.12]<-3

only works for two cases. If there are 3 cases (like in my case), the whole raster just fills up with the values from the last expression, in this case, 3.

What can I do to prevent this behaviour or is there any other method to replace the values and not have this effect?


texts: ArcScene 10.1 Orthophoto Visualization Issue

texts:
I am attempting to visualize an orthophoto in ArcScene 10.1 at the same resolution it is at in ArcDesktop 10.1. I have set the "Quality enhancement for raster images" to it's highest extent. But, the resolution is still low (i.e. it is very pixelated). In ArcDesktop I can see individual downed trees within the orthophoto, but within ArcScene I can barely see roads. Does anyone have any suggestions?

Thanks in advance.


texts:
I would like to use a basemap from the Ancient World Mapping Centre, available as maptiles via Mapbox, as basemap in my ArcMap desktop version 10.2.2.

It is possible to add a web tile layer to ArcGIS online as described in:

Using Stamen and MapBox tilesets as basemaps in ArcGIS.com

Is there a similar solution to add a custom map tile layer to ArcMap as in the question underneath?

Adding custom map tile layer using QGIS?

Note: I am looking for a solution that's not using external pluggin or ArcGIS Online.


texts: MapServer WFS service works in QGIS but not in ArcMap

texts:
I published a WFS service using MapServer. The source of the WFS layer is a FGDB.

When I add the service in QGIS, it works flawlessly.  When I add the service in ArcMap, all I can see is a layer named NoGeometry.

In ArcMap, I used 'Add Interoperability Connection', I chose the format WFS (Web Feature Service). In the Dataset, I put the the URL http://localhost:8080/cgi-bin/mapserv.exe?map=C:\ms4w\Apache\htdocs\aca\test.map

See the screenshot for the rest of the parameters.

When I drag the layer to add it to the map, I get this error: "This dataset appears to be invalid. Please verify that your settings or file contents are correct, and try again".

Where is the problem??


texts: Openlayers - Get extent of WFS vector layer and zoom to that extent when map

texts:
I would like to do the following:

1) While loading an OpenLayers map object get the lat/long Bounding Box of a WFS vector layer served from GeoServer.
2) Use the extent from the WFS vector layer to zoom so the entire layer is viewable.

This is the script I have now that uses the lat/longs from the WFS vector and then zooms. This script works fine.  I have a javascript function where the map and controls are created in the onLoad of the html Body tag. At the end of the script I call the loadZoomTo() function which then calls this line of script;

function zoomOnLoad() {
    this.map.setCenter(new OpenLayers.Bounds(-79.8961,41.1,-78.80809,42.2).getCenterLonLat(), 10);
}

In this script I have the lat/longs hard-coded but I need to be able to change these interactively depending on what option the user selects when he enters the application. If I remove the hard-coded lat/longs and try to replace them with something like this:

function zoomOnLoad() {
this.map.setCenter(new OpenLayers.Bounds(mLayers[a].getDataExtent()).getCenterLonLat(), 0);
}

it fails. I have also tried the event handler method as detailed here: OpenLayers, zoom to vector layer extent

but all I get is the html of the page with none of the OpenLayers stuff. No map, no controls etc..

Any suggestions would be appreciated.


texts: QGIS Masking part of objects at other layer

texts:
i'm using qgis 3.14. I have situation like in the screen below:
I'd like to hide that part of river which is under the bridge. So the result would be like:

I'd like to avoid editing original data set. Maybe it's possible by symbolization or masking, if not please tell the solution with editing data.


texts:
I had drawn a rectangle by adding element to graphics container like the code below:

IGraphicsContainer gc = mxDoc.FocusMap as IGraphicsContainer;
gc.Reset();
IRectangleElement rectangleElement = new RectangleElementClass();
IElement elem = rectangleElement as IElement;
gc.AddElement(elem, 0);

I want to remove the color from the rectangle or set the Transparency by 100, and keep the outline colored as it is.


texts: How to draw scatter plot from vector attribute table or vector (point) layer in GRASS?

texts:
I have a vector layer (point) with attribute columns cat, value1 and value2. I need to draw a scatter plot (value1 vs value2) in GRASS.

I can see scatter plot option for raster maps in GRASS? Is there a option available for vector/point map in GRASS?


texts: Dissolve without removing interior lines in ArcMap?

texts:
I'm fairly new to GIS and I'm having an issue that I hope someone can help me resolve?

Basically, what I'm trying to do is merge features based on a single attribute so that my table no longer displays the duplicates and the features are merged/grouped on screen. Other than going through my attribute table and manually selecting the features before clicking 'merge', the only other way that I could think to do this would be to use 'dissolve'. However, using this tool causes the features that are touching or overlapping to become one entity, which is the opposite to what I want to achieve.

In short, I wish to keep my features as separate entities that are grouped/merged together based on a single attribute whilst removing duplicates from my attribute table - is this possible?

Also, it's worth noting that I have 1000's of features to do this with, so I'm looking for an automated method of some kind?


texts: QGis not upgrading from 1.9.0. Ubuntu, Raring, using nghtly PPA

texts:
I'm currently running 1.9.0 on ubuntu raring, through the nightly ppas:

deb     http://qgis.org/debian-nightly raring main
deb-src http://qgis.org/debian-nightly raring main

But it doesn't update with software updater, even though I know there are more recent versions it should be pulling down. Can anyone think off-hand why that might be or what I could do to fix it?

I noticed that on the new qgis homepage, the nightly PPAs are not listed for raring... have they been discontinued for some reason, or temporarily broken?

Thanks

Simon


texts: Can I reference a way (street) in a node (representing a shop)?

texts: I was inspecting OSM wiki and some of dumps and I have a question now - suppose I want to add a node representing some shop to a particular street. Do I have to specify addr:street in that node? I see that that's what other users usually do. But that seems such of waste of effort, also prone to more human generated errors. Wouldn't it be more logical to simply reference an id of a particular way, representing that street in the given node?


texts: How much resources/performance is used by each additional projection available for a WMS service?

texts:
The service provider set up the GetCapabilities response for the WMS service. This defines which projections are supported by the WMS service and these are also the projections that can be requested.

I am wondering about the resources that are used when a large set of projections is made available.

I guess the map is stored with one projections and the GIS software on the server knows the formulas to re-project the data ?
But when is the image re-projected exactly in the process ?
And does this have any influence on performances and caching ?
What is the influence of the number of available projections on performances?


texts: I need to import LAS files to QGIS. I've never worked with GIS till now. Could anybody explain me this problem to the details? I've installed Quantum GIS v. 1.8.0.


texts: Aggregate items displaying only their number, and when zooming display each in QGIS

texts: Is there a way in QGIS to configure such a display like in https://velobs.2p2r.org/? For example, aggregated POI at low zoom level and zooming separates POI.


texts: "no current thread" bug when opening a file with gdal ogr (no multi threading)

texts:
This is a bug I cannot reproduce with simple code. I'm working with ubuntu 16.04 and gdal 2.2.3, in python3.

When I try to open a geojson with ogr in my programm:

src_ds = ogr.Open('/data/train/AOI_2_Vegas_Train/geojson/line/line_AOI_2_Vegas_img2389.geojson')

I get an error which is a segfault I guess:

Fatal Python error: PyThreadState_Get: no current thread

Current thread 0x00007f2bb0fcd700 (most recent call first):
  File "/home/myhome/project/env/local/lib/python3.5/site-packages/osgeo/ogr.py", line 7923 in Open

but when I try it in the python console, or in a small script, it works. I extracted the part of my code leading to this bug, but when I execute it separately from the rest of my code it works:

import gdal
import ogr

rgb_fn = '/data/train/AOI_2_Vegas_Train/RGB-PanSharpen/RGB-PanSharpen_AOI_2_Vegas_img2389.tif'
rgb_ds = gdal.Open(rgb_fn, gdal.GA_ReadOnly)
driver = gdal.GetDriverByName('GTiff')

def _get_raster_properties(raster_ds):
     # get reference properties
     x_size = raster_ds.RasterXSize
     y_size = raster_ds.RasterYSize
     geotransform = raster_ds.GetGeoTransform()
     projection = raster_ds.GetProjection()
     return x_size, y_size, geotransform, projection

x_size, y_size, geotransform, projection = _get_raster_properties(rgb_ds)

line_ds = driver.Create('/tmp/line_ds.tif', x_size, y_size, 1, gdal.GDT_Byte)
line_ds.SetGeoTransform(geotransform)
line_ds.SetProjection(projection)

line_fn = '/data/train/AOI_2_Vegas_Train/geojson/line/line_AOI_2_Vegas_img2389.geojson'
src_ds = ogr.Open(line_fn)  # this is failing in my code, but not here

gdal.RasterizeLayer(line_ds, [1], src_ds.GetLayer(), burn_values=[1])

this script is the equivalent of the code I am runnig in my program, with the same file treated.

Is this a gdal bug?


texts: Automatically backup from qgis to pgadmin

texts:
I have got postgres (pgadmin 3) + qgis link. And I think about next thing: how i can work with postgis layers, in a way that, backup file of my spatial database will created automatically for example every 10 minutes?
how can I set this up


texts: Deleting feature based on attribute using arcpy.da.Walk on file geodatabase

texts:
I'm very new to scripting and am having trouble figuring out how to use da.Walk to delete features based on attributes for my entire geodatabase. I got the script below to work on an individual feature class but haven't had any luck altering it to work on the whole geodatabase. I'm getting the following error:
RuntimeError: 'in_table' is not a table or a featureclass
 import arcpy
 arcpy.env.workspace = r'C:\Documents\Test.gdb'
 fc = arcpy.ListFeatureClasses()
 for files in arcpy.da.Walk(fc):
      for file in files:
           cursor = arcpy.da.UpdateCursor(fc, ["Zel"])
           for row in cursor:
                if row[0] != "AND":
                     cursor.deleteRow()


texts: OpenLayers 5.3 selected point feature does not return properties, but polygon does

texts:
I am having some issues in OpenLayers 5.3.0. The context is that I'm using Sharepoint lists as a backend to save feature data and OpenLayers is running in a Vue.js container. Everything works as expected:

The geoJSON feature strings are pulled in from SharePoint and aggregated into a feature collection.
As I load each feature I add custom properties - most importantly the ItemID from the SharePoint list.
The features load on the map.
I can draw new features and save all of them.

If at any point I get all features with code like

allFeatures = self.drawLayer.getSource().getFeatures();

the features will contain the additional properties I have added.

The problem I am having is that when use a select interaction the selected feature comes back differently if it is a point or if it is a polygon - using for example:

self.drawClickSelect.on('select', function(e) {
            self.selectedFeatures = e.selected;
            console.log(self.selectedFeatures);
          });

The feature comes back with the added properties if it is a polygon, but if it's a point the properties itemID and modified do not exist.

Examples
Polygon returns:

e {…}
dispatching_: (...)
disposed_: (...)
geometryChangeKey_: (...)
geometryName_: (...)
id_: (...)
listeners_: (...)
ol_lm: (...)
ol_uid: (...)
pendingRemovals_: (...)
revision_: (...)
styleFunction_: (...)
style_: (...)
**values_: Object
geometry: e
itemID: 13
modified: false**

Point returns:

0: e
dispatching_: (...)
disposed_: (...)
geometryChangeKey_: (...)
geometryName_: (...)
id_: (...)
listeners_: (...)
ol_lm: (...)
ol_uid: (...)
pendingRemovals_: (...)
revision_: (...)
styleFunction_: (...)
style_: (...)
**values_: Object
geometry: (...)**

self.drawLayer.getSource().getFeatures() returns: (Note that the properties I'm interested are present for each feature.)

"{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-9148122.009657301,4148932.8076870437],[-9148001.336280514,4149124.4626323017],[-9148174.732892036,4149232.323163778],[-9148225.35521765,4149153.4073993675],[-9148302.20055545,4149200.6003566994],[-9148335.00415841,4149155.742857091],[-9148551.727115743,4149293.291205446],[-9148561.789387284,4149278.5208011325],[-9148673.828665541,4149349.3877884443],[-9148735.371373663,4149252.6485635266],[-9148560.38595463,4149144.0397318634],[-9148567.23964962,4149133.2237431305],[-9148549.761260215,4149122.8749220027],[-9148544.360359373,4149132.5307491873],[-9148449.959080853,4149071.7602422866],[-9148456.755445173,4149061.5738832294],[-9148442.091697177,4149052.313675927],[-9148435.004098615,4149063.0075244764],[-9148340.775473515,4149000.2413342204],[-9148324.546920547,4149028.1325729615],[-9148199.787204698,4148944.980097804],[-9148184.701774493,4148972.1811943417],[-9148122.009657301,4148932.8076870437]]]},"properties":{"itemID":12,"modified":false}},{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-9149230.727863472,4149974.6889362033],[-9148917.825374642,4149776.7721247915],[-9149008.484043097,4149633.442597651],[-9149321.386531925,4149831.359409064],[-9149230.727863472,4149974.6889362033]]]},"properties":{"itemID":13,"modified":false}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-9149507.058888558,4150498.531335066]},"properties":{"itemID":21,"modified":false}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-9149383.916618763,4150420.6414931556]},"properties":{"itemID":22,"modified":false}}]}"

Is there some known issue with points that would make them different than polygons with regard to properties?


texts: Detecting checking/unchecking of rules in a rule-based rendering using PyQgis

texts:
Is it possible to programmatically detect checking/unchecking changes of a rule in a rule-based rendering ? I found no indication of existing signals ..

I got lost reading QgsVectorLayerRenderer Class Reference ...

Seeking preferably a solution for QGIS 3.x


texts: How can I see specific location on java desktop application?

texts:
I want to see a desired lat/lon point on the google maps on my java desktop application. I made some search so I think google map API is something which I looking for. Am I right? and Is there any useful example for me? (I am totally beginner at google Map API for java)
And I would not to use Android.


texts: Viewing OSM tag information of known element

texts:
I'm trying to filter out all the lakes in an area using osmfilter with natural=water and water=lake.  But some lakes are missing, so I want to view what tags are used for those missing lakes.
For example, I can see a lake graphic in the OSM map, but how can I view the tags of it?


texts:
I'm very new to QGIS, but have used Arc a bit. I've downloaded a shapefile of the Home Owners' Loan Corporation "residential security maps" from the 20th century (i.e. redlining maps). I collected the shapefile from the University of Richmond: https://dsl.richmond.edu/panorama/redlining/#loc=5/39.1/-94.58&text=downloads
In the QGIS browser pane I followed the filepath to the shapefile and dragged the holc_ad_data.shp onto the map pane. It will only show up as a single symbol (under "Layer Styling"). I want to show it (the layer's "HOLC_grade" field) as categorized or graduated symbols with the actual HOLC colors (red, yellow, blue, green), but whenever I change it away from a single symbol it disappears in the map pane.
Could anyone clue me in on how I can set the layer to the "HOLC_grade" field and change it from the "single symbol" style without having it disappear? I would also like to be able to have my mouse click on the polygons and bring up their attributes, instead of only being able to pan with the mouse.


texts:
I cannot download any of the plugins from the drop down menu, plugin installation failed. I have tried the re networking and also manual download but have not been able to find an answer for which folder to place ity into for this version of qgis. I am running Win10 pro.
Any help would be good.


texts: Filtering GeoServer GetFeatureInfo request

texts:
I have a 14-band raster stored in GeoServer. I use Jiffle based styles to calculate spatial indices on the fly and it works fine.
But when I clic on the map, GetFeatureInfo request returns the information of all the bands.
The request are like:
server/geoserver/wms/?REQUEST=GetFeatureInfo&SRS=EPSG%3A4326&BBOX=-56.12087982359703%2C-12.437183254414979%2C-56.12083982359703%2C-12.43714325441498&HEIGHT=594&WIDTH=528&LAYERS=layer&QUERY_LAYERS=layer&INFO_FORMAT=application%2Fjson&X=243&Y=263

And it returns:
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"id": "",
"geometry": null,
"properties": {
"Band1": 382,
"Band2": 343,
"Band3": 592,
"Band4": 321,
"Band5": 943,
"Band6": 4470,
"Band7": 5853,
"Band8": 5624,
"Band9": 6089,
"Band10": 6131,
"Band11": 2673,
"Band12": 1208,
"Band13": 3918,
"Band14": 0
}
}
],
"totalFeatures": "unknown",
"numberReturned": 1,
"timeStamp": "2021-03-24T11:18:15.298Z",
"crs": null
}

Is there any cql_filter or something that makes the operation and returns only the spatial index value?


texts: Combining two clips (feature classes) into one feature class in ArcGIS Desktop?

texts:
I clipped two shapefiles and created 2 clips but when adding the legend it adds all the same data in two clips but I want to add only one feature into the legend.

Is there any way to do it or can I merge two feature classes into one?


texts: QGIS: using postgis can I define a layer using sql?

texts:
I have a vector layer (call it master) that is stored in postgis and consists of a number of line segments.

I also have a collection of sub layers that consist of subsets of the master layer.  I can set up a 'link' table master_id, layer_name the says which segments belong to which sub layer -- some on the edges will be common.

Can I specify a layer using an join on the master and link tables?

Alternatively can I use a view at the postgres level to do this. I have never used views or stored procedures although I know the general principles.


texts: Getting 'module' object has no attribute 'ExtractValuesToPoints_sa' error in ArcPy while extracting values to points?

texts:
I am using ArcGIS 10.

I am trying to extract values to point using ArcPY but getting error as

  Runtime error : 'module' object has
  no attribute 'ExtractValuesToPoints_sa'

Below is what I am trying

import arcpy
from arcpy import env
from arcpy.sa import *

shape="G:\\Input\\ShapeFiles\\SA_Fishnet_0025_Points.shp"
raster="G:\\Input\\RasterFiles\\SA_Zonal_0025d.tif"
output="G:\\Input\\RasterFiles\\SA_Zonal_0025d_Zonal_PYTHON.shp"

arcpy.ExtractValuesToPoints_sa(shape, raster, output,"NONE","VALUE_ONLY")

I can do it without using ArcPy, but here I am learning it.


texts:
Is there a way to keep the outline of a building?

I am using a shp with lines showing the walls of a building (originally CAD file). I would like to keep only the outline of the building either as a line showing the perimeter of the building or as polygon showing the area of the building.

Any ideas?

EDIT

By outline I mean the boundary of the building
The example is just one floor of one building. I have multiple buildings per site and each has multiple floors. Each building_floor is a different map so there are hundreds per site.
The gaps is not a problem. Each gap is a door from a different layer. I can merge these two to close the gaps (basically close the doors).
I know I can manually draw a line or a polygon around each map but as I said it's a lot of work and I thought to ask first for a quicker way.
I already tried Concave Hull and Convex Hull but the results are not great - it needs fixing. Probably is better to start fixing this rather than start drawing each one from scratch.


texts: no popup menu with "right-click" on layer

texts:
I am trying to work with ArcGIS 10.4 on a MacBook Pro (OS X El Capitan 10.11.5) using Microsoft Windows via Parallels Desktop 11.

"Right-click" in this context normally = "control-click." No matter what combination I use, however, I can't seem to get the popup menu that should appear. I've also tried to find work-arounds for commands in that menu, e.g., "export data...," but can't figure out any functional workaround.

Does anyone have an idea how to solve the menu problem or know of an "export data..." workaround?


texts: Difference between ETRS89 / ETRS-TM32 and ETRS89 / UTM zone 32N

texts:
In QGIS I can choose among other things between "ETRS89 / ETRS-TM32" and "ETRS89 / UTM zone 32N". In my test project both projections displayed the layers at the very same place. So, where is the difference?

The proj4 codes are
+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs

What does +vunits=m mean?


texts: Buffer polygon inward by percentage of area

texts:
I would like to buffer a polygon inward by a percentage of its area, using the polygons below as an example.

How could I shrink each polygon inward resulting in a 35% reduction of area? I am using the term "reverse buffer" because I want the newly created polygons to have the same shape as the originals and be centered and within the originals, like the schematic below.

I am on ArcGIS 10.2 and have SA and ET Geotools.


texts: breaking road segments based on another shapefile in QGIS

texts: I am using QGIS 2.18 and have a road celterline and a grid shapefile which I am using to create an atlas book.  What I want to do is break each road at the edge of each grid that it crosses but I have been unable to find a way to do that.


texts: Can anyone suggest how I can change the length of an existing column (For ex. column name -- ABCD, column_Type -- Text, Length -- 254). It would be great to have solution on platforms like ArcGIS or QGIS.


texts: Setting up QGIS Server at low cost for small user community?

texts:
I want to setup QGIS Server at a low cost for a small user community (mainly me).

Parameters:

Some data is private so something like qgis cloud not an option (not enough money for paid option).
Primarily only me using this for experimentation. Sometimes I would like to share with a few other interested parties who may be scattered around the world.
Maybe one day I come up with something really interesting/useful and would like the option to scale up server capacity quite quickly in case it generates a lot of traffic.
Not too worried about 100% uptime - but equally don't want something that is falling over all the time. Essentially I can cope with the odd outage if it means its cheap.
Needs to be cheap.... ideally free.
Needs to be relatively simple to implement OR with a good tutorial. I have dabbled with this sort of thing before but still get a bit confused with too much server jargon & command line stuff - keen to learn at least!

QGIS server looks ideal because most of my work is in QGIS desktop and I want something that is relatively quick and easy to publish to.

Things I am considering:

Trying to setup QGIS Server on existing shared hostgator server... (but looking unlikely i can get sudo access).
Converting an old laptop into a home server.
Using Amazon EC2.
Try to setup QGIS Server on a bitnami stack and then use their cloud (maybe free option will work here?)
Finding a new low cost host option (could stretch to £10-20 a month at a push as a last resort)

Has anybody else faced a similar scenario? Can you suggest solutions that might work for me?

Many thanks in advance if so....


texts: Deriving bearing in degrees minutes seconds format from field calculator

texts:
I'd like to know if there's a formula for calculating bearings of lines in degrees minutes seconds format. I already have been able to get the bearings in decimal degrees using:

CONCAT(format_number(degrees(azimuth(start_point($geometry), end_point($geometry))), 2), '°')


texts: Can OGC Simple Feature polygon's interior boundaries share an edge?

texts:
Is the polygon pictured below compliant with the OGC Simple Feature standard?

There are two interior boundaries that share a segment so intuitively I would say it is not but I have not been able to find any confirmation in the documentation (in the 1.2.1 version, polygon is described on pages 26-28).

There are "rules that define a valid polygon" and I don't see there anything that would exclude polygons whose interior boundaries touch, but I also don't understand why, based on those rules, are some of the example polygons in Figure 12 invalid so I must be missing something.


texts:
We have SDE Running on SQL Server 2008. I have one user that receives this error once in a while, and not sure why he is the only one receiving it. I have not seen this error before. I checked permissions on the table it referenced (sdm_Equipment) and he has all the necessary permissions.

We have 4 people editing the database concurrently. We are all using ArcMap 10.3.1, and ArcGIS Server 10.3. The Error as seen below in the screen shot is:

Unable to save edits.
Underlying DBMS Error Cannot find the object "##SDE_Session3738_15 because it does not exist or you do not have permissions.][ CrescentLinkClarity.DBO.sdm_equipment]
The version has been redefined to reference a new database state.
[dbo.Default]

What would be causing this error? As I said above I checked the users permissions on sdm_equipment.


texts: Open access repository of general GIS spatial data?

texts:
As a newcomer to using GIS for research, I am wondering if there are online repositories for spatial data, similar to package repositories for Linux distributions.

For instance, I would be able to find things like "bathymetric map of the world's oceans", or "administrative regions map for the United States", etc.

Is there such a repository/website where I can download such data? Or do I have to "start from scratch" and manually draw a map of the world if I need one?


texts: How to add a sequence to the attribute table?

texts:
Is there a way to use the field calculator in QGIS to create a range of values in a count sequence with a specified input range?  For example, if I wanted to have the values 1 to 5 repeated within my table according to the first record in the table.

Thanks...


texts:
I have a table of Multiline Strings and need to generate a GeoJSON with the properties turning to EPSG 4326. I made the SQL below but is not in the correct format GeoJSON.

How do I create the GeoJSON using Postgis, so I have a FeatureCollection with the CRS 4326 set.


texts: GRASS GIS Segmentation issue on UAV image

texts:
I tried i.segment to perform segmentation on a part of UAV orthophoto image, but this always results in stripe-like rectangle segments which I do not expect. Could anyone figure out what happens?

Please let me note that r.info returns the size of this image as 6075 rows and 4525 columns, but g.region i -p returns only 86 rows and 7 columns.

The settings for image segmentation were as follows:

input: 3 bands (RGB) of raster maps
threshold: 0.01
minimum # of segment: 1
max. iteration: 100
others: default

3


texts: Projecting INEGI shapefiles from Mexico - not lining up

texts:
I have two INEGI shapefiles from Mexico (2000 and 2010) that I want to line up with one another. The 2010 file projects in the right location and I would like the 2000 file to match/line up with the 2010 one. The shapefiles contain various cities in Mexico - some line up almost identically (see example of Minatitlan) while others are dramatically different (see example below of Tetepango).

I originally posted a question about getting data from 2000 to project in the correct location in Mexico, which made the cities are least be in the generally correction location! See the response here: Projecting INEGI shapefiles from Mexico

What can I do to make the 2000 file match/line up with the 2010 one?

2010/2000 Data Source:
Projected Coordinate System:    Conica Conforme de Lambert
Projection: Lambert_Conformal_Conic
False_Easting:  2500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -102.00000000
Standard_Parallel_1:    17.50000000
Standard_Parallel_2:    29.50000000
Scale_Factor:   1.00000000
Latitude_Of_Origin: 12.00000000
Linear Unit:    Meter

Geographic Coordinate System:   ITRF92
Datum:  D_GRS_1980
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: need support on data driven pages in arcmap

texts:
I have searched for this answer online for 4-5 hours and am not finding the specific topic or solution to my problem.

I am creating a data driven pages mapbook of the state parks, state gamelands and state forests in Pennsylvania, with each page representing a county.  The problem that I am having is the locator (inset) map highlights a red border around the current extent of the main data frame in the locator map, rather than the individual county that the page is being created for.

I have clicked on the properties for the counties layer in the locator map's data frame, gone to the "Definition Query" tab and clicked on the "Page Definition" button, then checked the enabled checkbox, chosen "County_Nam" from the "Page Name Field:" dropdown and then clicked the "Match" bubble for "Show Features That," but this has not altered anything on my maps. The locator continues to highlight the extent, rather than the county.

Is there something that I am missing or another way to go about correcting this seemingly straightforward step in the data driven pages creation process?


texts: Geodatabase table as dynamic feature class domain

texts:
In ArcGIS Enterprise, is it possible to use a table stored in a geodatabase as a domain for a field in a feature class stored in the same geodatabase, such that it dynamically updates as new records are added to the table? Or do I have to use the table to domain tools each time I add new records to the table?

The idea would be to have a feature class of assets, and a table of owners, and restrict the values of the "asset_owner" field in the feature class to the entries in the table of owners


texts: Why is the Edit Vertices button (Draw Toolbar) enabled for some annotations but disabled for others?

texts:
I keep running into this:

Export Labels in Layer A & B to annotations (In Map, In current extent only)
Enable Focus Data frame
Click on annotation exported from Layer A
Click the Edit Vertices button - No problem
Click on annotation exported from Layer B
Edit Vertices button is DISABLED

What's the problem? When I export the same labels to geodatabase annotations, I can easily edit the vertices on any annotation.

Figure 1: Click on Layer A annotation - Edit Vertices button enabled

Figure 2: Click on Layer B annotation - Edit Vertices button DISABLED


texts: Problem transforming coordinates with PyProj

texts:
I would like to plot a marker using the Folium library but the data I have is in ETRS89 and Folium does not accept that projection. I am having problems transforming from EPSG25830 (ETRS89) to EPSG3857. I have read and followed step by step the tutorial in PyProj web and different ways of declaring the projections but something does not add up.
My code is:
from pyproj import CRS
from pyproj import Transformer

crs_in = CRS.from_epsg("25830")
crs_out = CRS.from_epsg("3857")

transformer = Transformer.from_crs(crs_in, crs_out)
print(transformer.transform(582040.535326,4.797240e+06))

(-221319.72875795467, 5361346.317915058)

That is weird because it should look like 43.259404, -2.026318 (not the same point, but very close) so I check the other way around:
  transformer_inv = Transformer.from_crs(crs_out, crs_in)
  print(transformer_inv.transform(43.259404, -2.026318))
  (834021.8587713321, -2.0147284128379286)

I do not understand what is happening or where I am wrong.


texts: Extract all raster cells outside a polygon

texts: How does one perform the very simple task of extracting all raster cell values outside a polygon? I've been looking through QGIS for hours with no solution to this simple problem.


texts: Use SearchCursor to populate SQL query from a table to run a selection on a feature class

texts:
I have a feature class with three fields of interest EV_CODE (text), TREE_SIZE_CODE (text) and CANOPY_COVER (number). I needed to find all possible combinations of those three attributes within the data set and then extract each possible combination as its own feature class.

I am creating a Python script that uses arcpy summary statistics to first create the table of all possible combinations (there are 158 possibilities in this data set). The values in the table all seem to be Unicode type (I don't know if this is also a problem, but it might be). Then I would like to use a SearchCursor with the table to create a SQL query for each row to select each possible combination from the feature class and export as its own feature class and a name that reflects the attributes selected.

The script works until it gets to the SQL query at the end where I keep getting ERROR 000358 that I have an invalid expression. I have included the error message below the script. Does anyone have a suggestion to fix the code so that the script will run?

# Find all possible combinations of three fields and select those areas to
# export as unique feature classes

# Import arcpy and set workspace
import arcpy
from arcpy import env
env.workspace = "C:/Users/astratton/Documents/Assignments/X_SideProject/TroutCreek/FSVeg_Combinations.gdb"
env.overwriteOutput = True

# Feature class
fc = "FSVegSpR06ExtractVegSIMPLIFIED"

# Identify fields of interest
field1 = "EV_CODE"
field2 = "TREE_SIZE_CODE"
field3 = "CANOPY_COVER"

# Use summary statistics cases to create table of unique cases
out_table = "Cases"
statistics_fields = [["OBJECTID", "COUNT"]]
case_field = field1, field2, field3
CasesTable = arcpy.Statistics_analysis(fc, out_table, statistics_fields, case_field)

# Add field delimiters to feature class fields
field1fc = arcpy.AddFieldDelimiters(fc, "EV_CODE")
field2fc = arcpy.AddFieldDelimiters(fc, "TREE_SIZE_CODE")
field3fc = arcpy.AddFieldDelimiters(fc, "CANOPY_COVER")

# Explore using cursor for Cases table to populate the SQL query
cursor = arcpy.da.SearchCursor(CasesTable, [field1, field2, field3])
for row in cursor:
    var1 = "'" + str(row[0]) + "'"
    var2 = "'" + str(row[1]) + "'"
    var3 = str(row[2])

    var1a = str(row[0])
    var2a = str(row[1])
    var3a = str(row[2])

    SQL = field1fc + " = " + var1 + " AND " + field2fc + " = " + var2 + " AND " + field3fc + " = " + var3
##    print SQL
    arcpy.Select_analysis(fc, "Unique_" + var1a + "_" + var2a + "_" + var3a, SQL)

Error Message...
Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.3\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\Users\astratton\Documents\Assignments\X_SideProject\TroutCreek\Python\FSVeg_Combinations_searchCursorTodayE.py", line 42, in <module>
    arcpy.Select_analysis(fc, "Unique_" + var1a + "_" + var2a + "_" + var3a, SQL)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\analysis.py", line 84, in Select
    raise e
ExecuteError: ERROR 000358: Invalid expression "EV_CODE" = 'None' AND "TREE_SIZE_CODE" = 'None' AND "CANOPY_COVER" = None
Failed to execute (Select).


texts: ogr2ogr - mapinfo to geojson. Merging shapes with certain values

texts:
I'm just starting out with this so apologies if I explain things poorly.

I have some mapinfo files (TAB, IND, MAP etc) and I am trying to convert to Geojson.  I can do this using ogr2ogr

ogr2ogr -f "GeoJSON" 1800-indoor.geojson coverage_LTE1800_Indoor.TAB

however the geojson is very large.  Actually however there is a lot of information in the mapinfo file (and thus the geojson file) that I don't actually require.  The random example tile of the map is below. Although it may not be apparent there are 5 'levels' of information ranging from 0 - 4.

When I convert the file all that is important is that I preserve the shapes that have a 0 value (although looking in the mapinfo file these values don't seem to be explicit, presumably they are encoded somehow?) but the shapes representing values 1,2,3, and 4 can all be merged together.  Visually therefore I would end up with two colours, one that represents 0 and another which represents 1,2,3 and 4 combined.

EDIT - more info added in response to comment:

list of files

example of data in data.map file

000a 0e15 0e16 1ba2 a6ab 1a2a 2fa5 a9b5
a7b5 d90f 1723 131f 2b0f 1723 4f57 639c
a4a9 a0a8 ada4 a8ad 161a 390d 1137 a5a9
b5a4 a8ad b2b6 dcbd bdf4 2b2f 55c8 ccd8
c7cb d0d3 d7fd c2c2 f900 0000 0000 0000
0000 0000 0000 0000 0000 0000 0000 0000


texts:
I am having an issue while trying to show the coordinates in the print layout in QGIS. While the coordinate format is in decimal and suffix, it is ok.

I need the format in degree, minutes, and suffix. But when I had selected that option, it was showing both north and south points in the latitude, and also the east and west in the longitude.

How to solve this problem?
There was a relative question like this one, but could not solve my problem.


texts: Creating a bottom-top aggregated slope model in QGIS

texts:
I have a Slope model and would like to create a bottom-top aggregated slope model across a waste dump. So I am looking at creating this through the slope's mean values.

One approach I thought of was to create a centerline of the waste dump polygon and use the SAGA cross profiles tool 1m apart. I would then create buffers 0.5m left and right of these lines (so each one lines up). These buffers could then be run through a zonal statistics process with the slope model to grab the mean slope value from bottom to top with the result a bottom-top aggregation around the waste dump.
The problem with this approach is the corners have gaps and overlaps and it is quite a poor alternative approach.

I need a solution like this:


texts: What does v.distance calculate excactly?

texts:
As I am from an economics department, I am quite a beginner in QGIS.
I read, that if I wanted to calculate the distance from a polygon to a point/line/another polygon, I should use the GRASS plugin v.distance.

Does anyone know how exactly the plugin calculates the distance? Let's say I want to know what's the minimum distance from Bavaria to the sea. Does v.distance take the center of Bavaria or does it divide the border into small points and then calculate the distance to the nearest sea spot?


texts: Custom QGIS linestyle expression function won't work with MultiLineStrings

texts:
After my experiences with a custom Python linestyle expression function (Wavy Edge on polygon in QGIS using geometry generators) I modified the code to support Polygons with holes, MultiPolygons and MultiLineStrings as well. See the following picture using the same Geometry generator styling for all layer types:

smooth(customLineStyle( $geometry,'45 30,45 0,45 -30,45 0'),5)

Unfortunately, QGIS crashes with MultiLineStrings in Windows 10 (QGIS 3.14.16-Pi).
I don't know, is it a bug or is there something wrong with my code (see below)?
import itertools, math, re
from qgis.core import qgsfunction,QgsPoint,QgsGeometry,QgsWkbTypes,QgsPointXY

def createCustomLine(geom,param):
    geom = QgsGeometry(geom)
    dst = geom.length()
    vertices = geom.asPolyline()
    start = (vertices[0].x(),vertices[0].y())
    end = (vertices[-1].x(),vertices[-1].y())
    l = [list(map(float,x.split(' '))) for x in param.split(',')]
    steps,offsets = zip(*l)
    d = sum(steps)
    newdst = dst / int(dst / d)
    newsteps = [x / d  * newdst for x in steps]
    d = sum(newsteps)
    cnt = int(dst / d)
    distances = list(itertools.chain.from_iterable(itertools.repeat(newsteps,cnt)))[0:-1]
    alloffsets = list(itertools.chain.from_iterable(itertools.repeat(offsets,cnt)))
    alloffsets.insert(0,0)
    distances = list(itertools.accumulate(distances))
    points2d = [(lambda g,d: (g.x(), g.y(),d))(geom.interpolate(d).asPoint(),d) for d in distances]
    distances.insert(0,0)
    points2d.insert(0,start) # prepend start point
    points = [QgsPoint(start[0],start[1])]
    for i,pt in enumerate(points2d[1:]):
        if distances[i+1] > distances[i]:
            corrAngle = -90
        else:
            corrAngle = 90
        qgsPt = QgsPoint(pt[0],pt[1])
        points.append(qgsPt.project(alloffsets[i+1],QgsPoint(points2d[i][0],points2d[i][1]).azimuth(qgsPt) + corrAngle))

    points.append(QgsPoint(end[0],end[1])) # append end point
    return QgsGeometry.fromPolyline(points)

def createCustomPolygon(geom,param):
    wkt = geom.asWkt()
    geometries = [QgsGeometry.fromWkt('LineString(%s)' % g) for g in re.findall('[^()]+',wkt) if g != ','][1:]
    geom = geometries[0]
    rings = geometries[1:]
    geom = createCustomLine(geom,param)
    geom = geom.convertToType(QgsWkbTypes.PolygonGeometry)
    if rings:
        for ring in rings:
            ring = createCustomLine(ring,param)
            ptArray = [QgsPointXY(pt.x(),pt.y()) for pt in list(ring.vertices())]
            geom.addRing(ptArray)
    return geom

@qgsfunction(args='auto', group='Custom', usesGeometry=False, referencedColumns=[])
def customLineStyle(geom,param,feature,parent):
    rings = None
    if geom.wkbType() == QgsWkbTypes.LineString:
        newGeom = createCustomLine(geom,param)
    elif geom.wkbType() == QgsWkbTypes.MultiLineString:
        wkt = "MultiLineString("+",".join([re.findall('^[^(]+(.+)',createCustomLine(part,param).asWkt())[0] for part in geom.parts()])+")"
        newGeom = QgsGeometry.fromWkt(wkt)
    elif geom.wkbType() == QgsWkbTypes.Polygon:
        newGeom = createCustomPolygon(geom,param)
    elif geom.wkbType() == QgsWkbTypes.MultiPolygon:
        wkt = "MultiPolygon("+",".join([re.findall('^[^(]+(.+)',createCustomPolygon(part,param).asWkt())[0] for part in geom.parts()])+")"
        newGeom = QgsGeometry.fromWkt(wkt)

    return newGeom


texts: What do I define the coordinate system of this data set to?

texts:
I have an excel spreadsheet with lat and lon columns, and 5 additional columns with associated data. I want to open the data points in ArcMap 10.2. I was told that the data was:

based on a Lambert conformal projection with the following specifications:

Central latitude/longitude 40N, 97W

First secant angle  33N
Second secant angle 45N

Origin at -1872000.0, 828000.0 from center.

How do I open this data in Arcmap with the correct spatial reference system?


texts: I selected my duplicates in my attribute table by Count(1,"ID")>1. Now I just only want to keep the duplicates with the highest areas. Is there a way to directly select them or a work around to keep them?


texts: Implementing WebAtlas_DE is failing due to Content-Type

texts:
I am trying to implement the German WebAtlas_DE.light WMS into one of our applications using the ArcGIS JS api.

The scenario is to load it in tandem when the dark-gray basemap is loaded. We have already implemented this with another WMS:

startup: function() {
    var self = this;
    this.dtk25 = new WMSLayer("http://sg.geodatenzentrum.de/wms_dtk250");
    this.dtk25.on("load", function(i) {
        var visibleLayers = self.dtk25.layerInfos.map(function(x) {
            return x.name;
        });
        self.dtk25.visibleLayers = visibleLayers;
    });
},

This is a snippet from our BaseMapManager widget, that handles all basemap configurations.
When we call the satellite base map, the dtk25 layer is loaded as well.
Trying the same with http://sg.geodatenzentrum.de/wmts_webatlasde.light however fails.

I've noticed that the Content-Type of the call to GetCapabilities is application/xml for dtk25 and application/vnd.ogc.wms_xml for the WebAtlas_de.light.

Using the Chrome Dev Tools I can replay that GetCapabilities call for http://sg.geodatenzentrum.de/wms_webatlasde.light?SERVICE=WMS&REQUEST=GetCapabilities and it downloads a file instead of opening it.

Printing the WMSLayer object to the console, after I've tried to load it, reveals a loadError-property:

"GetCapabilities request for http://sg.geodatenzentrum.de/wms_webatlasde.light failed. (Response is null.)"

What am I doing wrong here?


texts: Openlayers with allOverlays = true. How to create a dummy (base) layer?

texts:
I have an OpenLayers map with allOverlays set to true. Im adding 3 layers to the map. First, my own rendered version of OSM data (but only for levels 8-16), secondly an orthophoto layer (for levels 13-19) and 3rd layer is a vector, defined on all levels. Order is important, because of overlaying. Also, note that resolutions array i pass into the map object has no effect.

If i initialize the map in this way:

    map = new OpenLayers.Map('map_canvas',
    {
    projection: p3857,
    displayProjection: p4326,

    resolutions: [156543.0339, 78271.51695, 39135.75848, 19567.87924, 9783.939619, 4891.969809, 2445.984905, 1222.992452, 611.4962262, 305.7481131, 152.8740565, 76.43702827, 38.21851414, 19.10925707, 9.554628534, 4.777314267, 2.388657133, 1.194328567, 0.597164283, 0.298582142, 0.149291071, 0.074645535],
    units:'m',
    allOverlays: true,
    controls: [new OpenLayers.Control.Navigation(),
               new OpenLayers.Control.LayerSwitcher({'ascending':false})]
    });

...and add the layers in the order i specified, i have a problem, because of allOverlays: true it assumes the OSM layer as the fictional base layer (because its the first one added to the map) from which the map calculates its properties. These are incorrect because that layer isnt defined for all resolutions. I then simply copied the initialization code for my vector layer in front of XYZ layer (with OSM data) initialization and then it works ok, but the vector layer is rendered twice of course. I want to get rid of this dummy first layer which is needed to correctly initialize the map object. How can i do that?

I suppose i could use some world coastline shapefile and render it, but if its possible i'd rather just remove it somehow.


texts: Point data grouping based on location + huge dataset in QGIS

texts:
I need community help to solve my puzzle.
I have one huge dataset of points that represent places where a set of trucks stopped. Points go to the database if the reported speed of a truck was 0. The controlling system sends signal every couple minutes. It is ok for the driver to stop where others stop, but if there were just a couple of trucks stopped at a completely new place it is suspicious.

How would you select groups of points, which are not in the specific distance from any other group of points?  Result is selection of all points that have up to 10 (for example) neighbours but all these neighbours as a whole group have no other neighbours in criteria range.


texts: How to respond to marginal height error when using LiDAR DTMs in ArcGIS

texts:
I am comparing DTMs of a river valley over time using LiDAR data from the DEFRA download website. The Environment Agency state that this data has a relative height error of +/- 5cm.
The problem that I am having is that when I create a DTM of difference (i.e. minus one DTM from another) I find that the entire area has a net gain in elevation of approximately 0-3cm. Both DTMs are set to the same co-ordinate system and are alligned in the x-y dimesion. All of the other calculated elevation change seems appropriate (such as there is 1-4m elevation changes elsewhere). First of all, I am wondering if this is just the fact that the error margin is not precise enough to account for such small change? Unfortunately, it is making it look like the entire area has increased when in reality I don't think that it has. I realise 3cm is a tiny difference, particularly when the LiDAR has a resolution of 1m, but is there a way to fix this? I would prefer to have any values of change smaller than 5cm to instead measure as zero change.
Thanks in advance.


texts: How to use both Parse Path and Calculate Value in model?

texts:
Am I doing this wrong? I need to create hillshades for a set of raster, but I want the hillshades to retain the raster's original name. I used parse path for that, however there are names that are more than 13 characters so I used calculate value to edit the names. I used the expression r"%Value%".replace("-",""). It worked, I could see from the output that the value changed, however the hillshade does not adopt the naming even if I used %Value%.

Here is my model. I preconditioned everything so that the hillshade will be the last process.


texts: GeoServer Cluster - Malformed \uxxxx encoding

texts:
I’ve setup a Basic GeoServer cluster with a Load Balancer and a shared volume between the master and slaves with Docker. I’ve also got a Python script listening on changes to the data_dir (using watchdog), which then resets and reloads the slaves on changes on master.
This setup generally works great, but sometimes, from out of the blue, this error occurs for some layers claiming to be caused by some sort of malformed encoding:

Caused by: java.lang.RuntimeException: java.lang.IllegalArgumentException: Malformed \uxxxx encoding.
    at org.geowebcache.storage.blobstore.file.FileBlobStore.getLayerMetadata(FileBlobStore.java:691)
    at org.geowebcache.storage.blobstore.file.FileBlobStore.putLayerMetadata(FileBlobStore.java:647)
    at org.geowebcache.storage.blobstore.file.FileBlobStore.persistParameterMap(FileBlobStore.java:572)
    at org.geowebcache.storage.blobstore.file.FileBlobStore.writeFile(FileBlobStore.java:557)
    at org.geowebcache.storage.blobstore.file.FileBlobStore.put(FileBlobStore.java:477)
    at org.geowebcache.storage.CompositeBlobStore.lambda$put$6(CompositeBlobStore.java:182)
    at org.geowebcache.storage.CompositeBlobStore.lambda$readActionUnsafe$20(CompositeBlobStore.java:504)
    at org.geowebcache.storage.CompositeBlobStore.readFunctionUnsafe(CompositeBlobStore.java:486)
    at org.geowebcache.storage.CompositeBlobStore.readActionUnsafe(CompositeBlobStore.java:501)
    at org.geowebcache.storage.CompositeBlobStore.put(CompositeBlobStore.java:182)
    at org.geoserver.gwc.ConfigurableBlobStore.put(ConfigurableBlobStore.java:225)
    at org.geowebcache.storage.DefaultStorageBroker.put(DefaultStorageBroker.java:85)
    at org.geowebcache.layer.TileLayer.saveTiles(TileLayer.java:545)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getMetatilingReponse(GeoServerTileLayer.java:654)
    ... 111 more
Caused by: java.lang.IllegalArgumentException: Malformed \uxxxx encoding.
    at java.base/java.util.Properties.loadConvert(Properties.java:654)
    at java.base/java.util.Properties.load0(Properties.java:453)
    at java.base/java.util.Properties.load(Properties.java:404)
    at org.geowebcache.storage.blobstore.file.FileBlobStore.getLayerMetadata(FileBlobStore.java:689)
    ... 124 more

I can see that the metadata.properties-file for the GeoWebCache gets filled with sequences of \u0000, and the sequences seems to grow larger and larger, so I'm guessing the GeoServer finally runs out of memory, because after a while I can't request any tiles from the GeoServer.
I've only found one other post about this problem here: https://osgeo-org.atlassian.net/browse/GEOS-9272
Does any one know why this problem occurs and how one could fix it? I'm using GeoServer 2.17


texts:
I´m stuck with getting GWC working for one layer. Geoserver version is 2.9.
The layer in question has

Native Coordinatesystem: EPSG:4326
MinX: 19.0583763  MinY: 43.1548072  MaxX: 19.165252  MaxY: 43.2070563

I´ve created a GWC gridset based on EPSG:4326 and calculated the gridset borders from the scope of the reference system:

MinX: -180  MinY: -90  MaxX: 180  MaxY: 90

I´ve added the gridset to my layer. On saving I already get warnings the the bounds do not overlap with a gridset. When visiting a tile and checking the return header I see:

geowebcache-cache-result:MISS
geowebcache-miss-reason:

and geoserver logs show:

2018-03-09 12:34:38,493 INFO [org.geoserver.gwc] - Error dispatching tile request to GeoServer
org.geowebcache.grid.OutsideCoverageException: Coverage [minx,miny,maxx,maxy] is [128, 128, 127, 127, 8], index [x,y,z] is [141, 162, 8]

What is going wrong?


texts: Defining projection in arc toolbox (JSON to shapefile) in python

texts:
I have a arc toolbox which converts JSON txt file to a shapefile but we have to define the projection.

How can I hard code the defined projection so that user don't have to input the projection ?

I need to define projection parameter to web mercator and to add a function to reproject the shapefile to NAD 83 UTM.

I am editing in Python.


texts: Iterate network analyst closest facility, arcpy error

texts:
I am working on simulating resource deployments within arcpy network analyst.
I have individual incidents stored in a file database (they have to remain single because later on the script must do a check to ensure there is a resource to deploy. My problem is with the iteration. The first closest facility analysis works but then the script fails before completing the next one.
It returns

arcgisscripting.ExecuteError: A python error occurred.
Failed to execute (FindClosestFacilities).

This is my code snippet below

import arcpy
import time
arcpy.env.workspace = r"c:\Simulation\incidents.gdb"
arcpy.env.overwriteOutput = True
start = time.clock()
arcpy.CheckOutExtension("Network")
featureclasses = arcpy.ListFeatureClasses()

print featureclasses

for fc in featureclasses:
    arcpy.env.workspace = r"c:\Simulation\incidents.gdb"
    arcpy.env.overwriteOutput = True
    inNetworkDataset = r"c:\Simulation\tempFGDB.gdb\deccan_street\deccan_street_ND"
    inIncidents = fc
    inFacilities = r"c:|Simulation\proposed_csl.shp"
    outGeodatabase = r"c:\Simulation\NA_outputs.gdb"
    outRoutes = "Routes"
    outDirections ="Directions"
    outClosestFacilities = "ClosestFacilities"
    measurement_units = "Minutes"
    with arcpy.da.SearchCursor(fc,("DATE_TIME", "OBJECTID")) as cursor:
        for row in cursor:
            vALUE = str(row[1])

    # Run FindClosestFacilities. Choose to find only the closest facility.
    arcpy.na.FindClosestFacilities(fc, inFacilities, measurement_units, inNetworkDataset, outGeodatabase, outRoutes + "_" + vALUE, outDirections + "_" + vALUE, outClosestFacilities + "_" + vALUE, Number_of_Facilities_to_Find=1)

print "Script completed successfully"

elapsed = (time.clock() - start)
print("Execution time: " + str(elapsed))


texts: MULTIPOLYGONM vs MULTIPOLYGON with shp2pgsql

texts:
I have already seen this answer

I am running the following:

shp2pgsql -I -s 3401 -t 2D theShapeFile.shp table_name > file_name.sql

And the response I see is:

Shapefile type: PolygonM
Postgis type: MULTIPOLYGONM[2]

My database does not support MULTIPOLYGONM (and I cannot upgrade).

How can I force the output to be MULTIPOLYGON instead of MULTIPOLYGONM ?


texts: How to set/update columns by using a trigger?

texts:
I want to set/update attributes (alignment_horizontal, alignment_vertical) of a PostGIS geometry (point) by using a trigger.

CREATE TABLE point (
gid serial NOT NULL,
geom geometry(point, SRID),
label varchar(255),
x_cord double precision,
y_cord double precision,
alignment_horizontal varchar(255),
alignment_vertical varchar(255),
CONSTRAINT point_pkey PRIMARY KEY (gid)
);

CREATE TRIGGER label_alignment_trigger BEFORE INSERT OR UPDATE
    ON point FOR EACH ROW
    EXECUTE PROCEDURE label_alignment_triggerfn();

Is there any way to to set/update two columns in one trigger function?

CREATE OR REPLACE FUNCTION label_alignment_triggerfn()
RETURNS trigger AS $body$
    BEGIN
    NEW.alignment_horizontal :=
        CASE
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 22.5 THEN 'Center'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 112.5 THEN 'Left'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 202.5 THEN 'Center'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 292.5 THEN 'Right'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID)))  >=  292.5 THEN 'Center'
        END;
    NEW.alignment_vertical :=
    CASE
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 22.5 THEN 'Bottom'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 112.5 THEN 'Half'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 202.5 THEN 'Top'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID))) < 292.5 THEN 'Half'
            WHEN  degrees(ST_Azimuth(NEW.geom, ST_SetSRID(St_MakePoint(NEW.x_cord, NEW.y_cord), SRID)))  >=  292.5 THEN 'Bottom'
        END;
    RETURN NEW;
    END;
$body$ LANGUAGE plpgsql;

Edit: Trigger function edited.


texts: Export only select feature properties to a CSV file from Earth Engine

texts:
I want to calculate the number of fires for a region and export the result to Google Drive as a CSV file. The code I'm using to do this seems to work, but the object fire_ind_count contains many additional properties besides fire count. How can I alter my code to obtain an object with only the count of fires to export?

// Load country shapefile
var lsib = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
var vn_shape = lsib.filterMetadata('country_na', 'equals', 'Argentina');
print (vn_shape);

// Load fire counts image
var fire = ee.ImageCollection('FIRMS')
             .filterBounds(vn_shape)
             .filterDate('2019', '2020');
var scale = fire.first().projection().nominalScale();

// Filter fire with more than 50% confidence and add a new band representing areas where confidence of fire > 50%
var filterConfidence = function(image) {
  var line_number = image.select('line_number');
  var confidence = image.select('confidence');
  var conf_50 = confidence.gt(50).rename('confidence_50');
  var count_band = line_number.updateMask(conf_50).rename('count');
  return image.addBands(count_band);
};
var fire_conf = fire.map(filterConfidence);
print('fire_conf', fire_conf);

// Count for individual image.
var countIndividualImg = function(image) {
  var countObject = image.reduceRegion({
    reducer: ee.Reducer.countDistinct(),
    scale: scale,
    geometry: vn_shape,
    maxPixels: 1e9
  });
  return image.set(countObject);
};
var fire_ind_count = fire_conf.map(countIndividualImg);
print('fire_ind_count', fire_ind_count);
Export.table.toDrive({
  collection: fire_ind_count,
  description: 'FireCountsInAustralia_gt50_conf'});


texts: Adding point data from text file when coordinates are in OSGB36

texts:
I have a text file in which coordinates are in the OSGB36 National Grid reference system used in Great Britain, so the coordinates are not in longitude and latitude (e.g. coordinates look like this: SD481731, and are only in one column).
I would like to display these points in QGIS 2.4.0.
However, when adding a delimited text layer, it is mandatory to specify the X and the Y fields.

Can I circumvent this and display OSGB36 coordinates imported directly from one column in a text file?


texts: How do I classify a highly heterogeneous built-up area?

texts:
I am using Landsat 8 imagery for supervised classification. I have used maximum likelihood classifier to classify separate classes such as forest, grassland, agriculture, plantation, water body and built-up area. The built-up areas in my study site are small towns with a lot of green cover, and therefore, had been classified as either one of the vegetation classes stated above.

Should I create a polygon and clip those areas where these towns are located and just classify them as built-up area? Also, please note that I will be doing change detection on this classified image and don't know whether it would be advisable to use this method to 'manually' classify these built-up areas.


texts: gdal_merge - creating intermittent (null seams) in mosaics

texts:
I’m supporting the web services development at Geoscience Australia's Digital Earth Australia (DEA) project

The DEA landsat coverage for Australia is stored as  > 1400  1 degree tiles in an Australian Albers projection.

Mosaicked together a 3 band byte file is > 112 GB.  Its a big data set and there are multiple continental products to be published.

We’ve been using the NCI to process these data – and have run into problems creating national mosaics and pyramids for our web services products.

We originally found that running ‘retile’ on the full national mosaic was taking > 14 hours on the NCI - and it becomes very unstable when rendering the higher resolution pyramids and crashes more often than not.

So as an alternative workflow .. we moved to a processing model based 16 degree processing blocks across the continent.

Using these tile groups we use ;

Gdal_buildvrt (used to create a virtual mosaic)

     *gdal_buildvrt(l_vrPathName,l_targetDir):
    cmd = r'gdalbuildvrt ' + l_vrPathName + ' ' + l_targetDir + '/*.tif'**

Gdal_translatefiles (create a byte precision mosaic)

    *gdal_translateFiles(l_vrPathName,l_inputEPSG,l_tempPathName):
    cmd = r'gdal_translate --config GDAL_CACHEMAX 500 -a_srs "' +     l_inputEPSG + '" -scale 0 4095 0 255 -b 1 -b 2 -b 3 ' + '-of GTiff ' + '-ot Byte ' + '-co "TILED=YES" ' + '-co "BLOCKXSIZE=512" ' + '-co "BLOCKYSIZE=512" ' + l_vrPathName + ' ' + l_tempPathName*

Gdal_warp (reprojects mosaic to wgs84 geographic’s)

   *gdal_warp(l_tempPathName,l_inputEPSG,l_destinationPathName):
   cmd = r'gdalwarp -r ' + g_samplingMethod + ' -srcnodata ' + g_srcNoData + '-of GTiff -multi --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS=ALL_CPUS -s_srs ' + l_inputEPSG + ' -t_srs "' + g_outputEPSG + '" ' + l_tempPathName + ' ' + l_destinationPathName*

Gdal_retile ( creates 7 layer pyramid structure)

      *gdal_retile(l_destinationDir, l_destinationPathName ):
      cmd = r'gdal_retile.py -v -r near -ps 2048 2048 -levels 7 -pyramidOnly -co "TILED=YES" -co "COMPRESS=LZW" -targetDir ' + l_destinationDir + ' ' + l_destinationPathName*

Gdal_merge (creates mosaics from pyramid structure)

    *gdal_merge(mosfile, in_list):
   cmd = r'gdal_merge.py -n 0 -v -o ' + mosfile + ' --optfile ' + ' ' + in_list*

This model proved quite fast as it lent itself to parallel processing and is very stable.

BUT what we are found is that ‘merge’ is creating small intermittent seams (no data) between the tile blocks.

For testing I ran a workflow with a full 1 degree overlap between the tile blocks .. but the seams (albeit more intricate) still persist.

I suspect the gdal_warp is the initial source of the problem .. but the coverage is so large using the australian albers projection in geoserver just doesn't work.

Try as I might - I can't find a gdalwarp command line structure using; -tr -tap  that works!?

gdal_merge (in its current form),  works only on the neat lines of the input tiles and ignores any areas of overlap.

Question 1 : Is there a gdal mosaicing pathway that supports 'feathering' or a NN mapping in areas of overlap ?

I'd welcome any comments or suggestions on an alternative gdal workflow.


texts:
I have a catalog of just over 3000 LAZ files that I want to split into regions using clip_rectangle().  The total catalog size is 60GB.  My PC has 64GB of RAM but I get an error that lidR is unable to create a 2.8GB vector when I try to create my subset.
Is the simple answer to get more RAM or is there a way to avoid this without having to add more RAM?
Here's my script:
library(lidR)
set_lidr_threads(6L)
ctg <- readLAScatalog("c:/noaa/utm/ground/")
#Make sure we plot the catalog so we can see where we want our clip regions
plot(ctg)
subset = clip_rectangle(ctg, 5550000, 2795000,570000,27972500)


texts: MapServer map file shows Raster image as white

texts:
I'm using PAN image for raster publishing but the output in browser shows a full white image? Please give me the solution?
My map file is ass follows:

MAP
NAME "terra"
# Map image size
SIZE 600 600
UNITS dd
EXTENT 77.0962025047 23.8718740333 84.6389301196 30.4075924299
  PROJECTION
   'proj=longlat'
   'ellps=WGS84'
   'datum=WGS84'
   'no_defs'
  END
  # Background color for the map canvas -- change as desired
  IMAGECOLOR 155 155 155
  IMAGEQUALITY 95
  IMAGETYPE png

  OUTPUTFORMAT
   NAME png
   DRIVER 'GD/PNG'
   MIMETYPE 'image/png'
   IMAGEMODE RGBA
   EXTENSION 'png'
  END
  # Legend
  LEGEND
      IMAGECOLOR 255 255 255
    STATUS ON
    KEYSIZE 18 12
    LABEL
      TYPE BITMAP
      SIZE MEDIUM
      COLOR 0 0 0
    END
  END

  # Web interface definition. Only the template parameter
  # is required to display a map. See MapServer documentation
  WEB
    # Set IMAGEPATH to the path where MapServer should
    # write its output.
    IMAGEPATH "../tmp"

    # Set IMAGEURL to the url that points to IMAGEPATH
    # as defined in your web server configuration
    IMAGEPATH "../tmp"

    # WMS server settings
    METADATA
      'wms_title'           'terra'
      'wms_onlineresource'  'http://my.host.com/cgi-bin/mapserv?map=wms.map&'
      'wms_srs'             'EPSG:4326'
    END

    #Scale range at which web interface will operate
    # Template and header/footer settings
    # Only the template parameter is required to display a map. See MapServer documentation
  END
  LAYER
    NAME "UPC1PAN10SISDP2012"
    TYPE RASTER
    OFFSITE 0 0 0
    DATA 'Z:\Rasters\10K\10D\UP\UPC1PAN10SISDP2012.img '
    METADATA
      'wms_title' 'UPC1PAN10SISDP2012'
      'WMS_ENABLE_REQUEST' '*'
    END
    STATUS DEFAULT
    TRANSPARENCY 100
    PROJECTION
    'proj=longlat'
    'ellps=WGS84'
    'datum=WGS84'
    'no_defs'
 ''
    END
  END

END


texts: Problem with CORS on GeoServer running on Apache Tomcat

texts:
I am currently working on a little GIS web app using MapBox on the frontend and GeoServer as my WMS. I have GeoServer 2.15.0 running under Apache Tomcat 9.0.20. Both of them run on localhost but on different ports.

Problem is, when I want to get some data from GeoServer to load it into my app I get a CORS response.

I have followed many threads and the Apache Docs on how to enable CORS both for Tomcat and for GeoServer using the XML filter configuration. I've followed many answers but nothing seems to work.

As a side note, this only happens to me on Windows, at first I tested it on MacOS on my office computer and everything worked just fine on the first try. So then when I got home and repeated the process on Windows, it doesn't work, and I don't see what else I am missing.

On MacOS i installed Tomcat using Homebrew, on Windows using the Windows Service Installer. I installed GeoServer using the .war file.

This is the error I get on the console and the headers of the request. I am not sure if this is the best way to get JSON data from GeoServer, but it worked fine on the Mac. Following the recommendation, I deactivated Windows Firewall to test it but it still fails.


texts:
I am using QGIS 2.14.3.
I have a Layer with different attribute string columns. I've colored one layer for the first column with properties->style and random colours (about 900 different colours). I duplicated the layer and colored it for the second column.

In the first and second column are the same names, but they are differently distributed and not all rows in the second column are filled. Now I want to color the second column so that the same names have the same colour as in the first column.

Is there an easy way to do that?


texts: How to export OpenStreetMap data from QGIS for printing

texts:
For my work I want to print a road map of Mallorca on A0 size. Therefore it is important that the base file is a vector file. After doing some research I found that I should be able to open OpenStreetMap data in QGIS and from there to export a CSF or PDF file. But so far I have only just been able to open the OpenStreetMap in QGIS.

Which plugin should I use to select only part of the osm data and how can I export the selected part?


texts: Moving and converting virtual fields in QGIS

texts:
I have two questions about virtual fields:

Is it possible to move a virtual field? I'm using Table manager to move "normal" fields but virtual fields doesn't show.

Is there any simple way to convert a virtual field to a normal field? Or even converting all virtual field in a layer at once?


texts: Create DEM from either few points or contours without elevation attributes

texts:
I have DWG file that contains couple of layers, among others I have "elevations" layer (that consists of only 80 points) and "contours" layer that can be seen on the images below. My goal is to make DEM from any resource that would be suitable to accomplish this task.

In the above image you can see Elevations + borders layer, and in the below one can see Elevations + contours + borders layer.

I don't have any CAD installed so I am using online software to read file and also to convert it to dxf. When I converted file to dxf and opened it in QGIS I could only see elevations of points but not contours in attribute table (I assume that even this bigger black elevation on the image above 91.00 is missing). Any kind of interpolation that I was doing in QGIS (How can I create a DEM from Point Data) didn't gave me result which would be similar to the contours that I already have in dwg file. I was also trying to do something in Surfer v.8 at first and I wasnt satisfied with result at all, and in Surfer v.13 I was given a notice that I have inssuficients number of points. I guess the real problem is that I have really only 80 points that dont represent good enough terrain but maybe I am wrong. Then I am wondering how did the dwg file creators made contour lines.

Should I continue searching for solution to this problem or should I start searching how to join elevation attributes to the contour layers (by calculating them manually or any other approach that anybody suggest) and try to make DEM from it?


texts: ArcGIS Javascript API - switch case to display attributes in a popup

texts:
var vliTemplate = new InfoTemplate("Visual", "
VLI Number: ${VLI_VLI_NO}
	"+

"

Not Greened Up CFLB:${Tbl_ResultSummary_VQO_Not_Greened_Up_Ha:NumberFormat(places:0)} Ha
	"+
"

Greened Up CFLB:${greened:compare} Ha
	"+
"

Maximum Target Percent:${Tbl_ResultSummary_VQO_Pct_tgt:NumberFormat(places:0)}%
	"+
"

Maximum Target:${Tbl_ResultSummary_VQO_Target_Ha:NumberFormat(places:0)} Ha
	")

compare = function (key) {
  var result = diff;

  switch (key) {
    case "greened":
      diff = data.Tbl_ResultSummary_VQO_HA - data.Tbl_ResultsSummary_VQO_Not_Greened_Up_Ha;
      break;
  }

  return result;
};

I've been playing around with the code in my switch case function trying to get wording or placements of words to work for me, but I can't seem to get my popup to display anymore. In the popup, I want the line "Greened up CFLB" to show the calculation of Tbl_ResultSummary_VQO_HA - Tbl_ResultsSummary_VQO_Not_Greened_Up_Ha but now the popup won't even display (because I know I'm not writing my fucntion/switch statement correctly)

Can someone please provide feedback on where I'm going wrong in my code? I know its a simple case of not having the right things in the right place...but I'm at a loss.

Thanks!


texts: Invalid variable for Set LAS Class Codes Using Features in ModelBuilder of ArcGIS Pro

texts:
I am attempting the following process in Model Builder in ArcGIS Pro 2.5. I have an original .las file as input into Extract LAS. What I failed to mention initially is this is NOT LIDAR, this is LAS from UAS Photogrammetry. While I have completed the entire process I am about to describe without issue, manually, I am now having a lot of difficulty stringing these tools together.
When run alone, Extract LAS works, but I cannot work out how to string it into Set LAS Class Codes Using Features.  Even using 'Out' as a precondition, running the tool yields 'Parameters are not valid.  The value does not exist'.
From what I know, Extract LAS generates a new .las and a .lasd, however only the .lasd displays as output in Model Builder.  This is problematic for my next tool, because it seems Set LAS Class Codes Using Features only accepts a .las as input.  For example, it will accept my original .las file as input, but of course I want the output version that is generated from Extract LAS.
I have tried using Path and Name elements from Parse Path as input directly into Set Las Class Codes, I have tried creating a new variable from Path and Name to use as input, I have tried Calculate Values.  For each, I've tried every combination of stringing these elements together:  Name alone, Name with the .las extension, Path and Name, Path and Name with a .las extension, etc etc. but all are invalid input.  I have tried using 'Out' as a precondition.
Can anyone see what I need to do to string this together?
I've tried too many combinations to list here but here is one example.  Currently, the original .las input is turned on in my map.

After posting the above, I ran Extract LAS, then added the Set LAS Class Codes Using Features Tool and navigated directly to the new Extracted .las file I needed as input.  When I re-ran the model, it set the Input LASD Dataset as a Precondition (shown at bottom) and this worked.  I am positive when I tried this, on my own, the model still threw the error indicating the vale did not exist.
Is there some trick to setting preconditions that I might be missing?
Even more, now that I have a Reclassified LAS Dataset, the next tool now will not let me connect to it as the next input, nor will the LAS to Multipoint tool let me input any of my LAS files even in a new model when LAS To Multipoint is the only tool.  The only options for input seem to be the Output Folders from the first step, or to set the Reclassified Dataset as an Extent or Precondition when I try to connect it with an arrow.  I can only add the LAS as input if I manually navigate to a folder and select it.
Is there no way to input my individual LAS file as input from within ModelBuilder, or must I put the Reclassified LAS output into it's own folder in order to connect it to LAS to Multipoint?


texts: Can't add a specific WMS on Leaflet which worked on QGIS

texts:
I am learning to use Leaflet. I tried to add a layer from  https://maps.nlsc.gov.tw/OpenData/wms. My code looks like this:

var wms = L.tileLayer.wms("https://wms.nlsc.gov.tw/wms", {
    layers: 'EMAP2',
    format: 'image/png',
}).addTo(map2);

But the layer didn't show on the map. When I checked what is being returned by the developer's tool, I found they were xml files and not images! So this may be the problem, but I don't know how to fix it.
Other things I've tried:

Checking with another layer, "CITY". Doesn't work. Still no layer is showing.
Check if something is really wrong with my code by trying to load a layer from another source by changing the url to "http://ows.terrestris.de/osm/service", and layer to "OSM-WMS". This works!
Check if something is something is wrong with the coordinates by loading the 'EMAP2' from the first source and "OSM-WMS" layer from the second source in QGIS. I set the CRS to EPSG:3857 for both. This works! They overlay just fine.
Check what request the browser is sending using developers tool. For the "EMAP2" layer, this is a sample of what is being requested (as there are multiple request):https://wms.nlsc.gov.tw/wms?service=WMS&request=GetMap&layers=EMAP2&styles=&format=image%2Fpng&transparent=false&version=1.1.1&width=256&height=256&srs=EPSG%3A3857&bbox=13511620.615914037,2886262.188048257,13521404.55553454,2896046.1276687565


texts: Saving each file in for loop gives ERROR 000872?

texts:
I am beginner to Python programming. I know, I am asking a really silly question. I am trying to extract DEM for each watershed from a bigger DEM. Therefore, I am using extract by mask. The bigger DEM is in raster format whereas the watershed boundaries are polygon. I have the following python code which works well. However, I am unable to save each watershed DEM after extraction in the loop. using my code I can only save the 1st watershed,then it says"ERROR 000872: Output Raster: Output C:\Subhasis\Project-Biological_monitoring\LIDAR-DEM\i exists. It cannot be overwritten since overwrite is off".

Code:

import arcpy, os
from arcpy import env
from arcpy.sa import *

#set work environment
env.workspace = "C:/Subhasis/Project-Biological_monitoring/Individual_watershed/Watershed_Final"

outws="C:/Subhasis/Project-Biological_monitoring/LIDAR-DEM"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

#List shapefile
fcList = arcpy.ListFeatureClasses()

for i in fcList:
    inRaster="nj10ftffil"
    outExtractByMask = ExtractByMask(inRaster, i)
    #outname=os.path.join(outws,i)
    #outExtractByMask.save(outname)
    outExtractByMask.save("C:/Subhasis/Project-Biological_monitoring/LIDAR-DEM/i")


texts: GRASS init error when adding GRASS as Plugin in QGIS 3.18.2

texts:
When I add GRASS from a list of Plugins, I receive an error:

GRASS init error: GRASS was not found in 'C:/PROGRA~1/QGIS
3.18/apps/grass/grass78' (GISBASE), provider and plugin will not work.

Although the Plugin has been installed, when opening the GRASS drop-down menu, Open mapset etc. is greyed out.
I received the same error after installing the stable version 3.16.6.
I have so far tried to uninstall and reinstall various versions and I'm unable to get this going. Windows Version 10.0.19042 Build 19042.


texts:
A standard PostGIS point-in-polygon query might look like this. (example tutorial)

SELECT p.*,b.polyname
 FROM points p, shapes s
 WHERE ST_Within(p.the_geom,s.the_geom)

Imagine this was not a spatial join. Then it would look like this:

SELECT i.*,c.categoryname
 FROM items i, categories c
 WHERE i.zipcode = c.zipcode

Most SQL folks would advise rewriting that query using a variety of the JOIN keyword syntax:

SELECT i.*,c.categoryname
 FROM items i
 JOIN categories c ON i.zipcode = c.zipcode

Does PostGIS offer this sort of syntax for spatial joins? I can't find it. All I see are a lot of functions that return true or false, leading to the workaround syntax. Ideally, I would think a spatial join would look like:

SELECT p.*,b.polyname
 FROM points p
 JOIN shapes s
 ON
 p.the_geom IS WITHIN s.the_geom


texts: How to do raster conditional comparisons in .NET?

texts:
I am comparing a raster with a constant value by raster calculator. Now I want to do this in a c# code.

Question

Is there any process to do this kind of raster conditional comparison in .net?


texts: Why segmentation fault when QgsVectorlayer not deleted?

texts:
Config: Ubuntu 20.04, QGIS 3.10, Python 3.8.5
My pyqgis scripts are doing the job, but some of them end with a segmentation fault when executing the closing instruction .exitQgis().
The segfault appear each time that I instantiate a QgsVectorLayer that does not get deleted before the end of the script.
Adding a del my_vectorlayer before the end solves the segfault. In the following example it is a bit trivial, but on longer processes using classes, such detail can be easily forgotten as I'm not used to delete explicitly variables.
Is it normal?
Is there another way of doing it better?
I have been looking in documentation but didn't find a proper explanation.
# -*- coding: utf-8 -*-
import sys
sys.path.insert(0, '/usr/lib/python3/dist-packages')

from qgis.core import (
    QgsApplication,
    QgsVectorLayer
)

print('start app')
QgsApplication.setPrefixPath("/usr", True)
myapp = QgsApplication([], False)
myapp.initQgis()

test_vector_file = 'test.gpkg'

print('process layer')

my_layer = QgsVectorLayer(test_vector_file, 'test layer', 'ogr')

# if del instruction omitted => Segmentation fault (core dumped)
del my_layer

myapp.exitQgis()
print('closed')

Segmentation fault (core dumped) happens after the print('closed')