texts: R - Spatial Join Between SpatialPoints (GPS coordinates) and SpatialLinesDataFrame

texts:
I am working on a university project which combines data science and GIS. We need to find an open-source solution capable of obtaining additional information from a massive GPS coordinates dataset. Clearly, I cannot use any API with daily request limit.

Here you can find a sample of the dataset the Professor provided us:

longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)
ID <- seq.int(1, 10)

The first step was joining my SpatialPoints with a SpatialPolygonsDataFrame using over() of rgeos. The SpatialPolygonsDataFrame was obtained through getData('GADM', country='ITA', level=3) of rgeos.
For this first accomplished task, the objective was to associate to each GPS coordinates the information about City and Region which they belong to.
An example of the result I was able to obtain is:

require(sp)
require(rgeos)
my_spdf <- SpatialPointsDataFrame(coords = longlat, data = ID, proj4string = CRS(" +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 "))
italy_administrative_boundaries_level3 <- getData('GADM', country='ITA', level=3)
result <- over(my_spdf, italy_administrative_boundaries_level3)[, c("NAME_0", "NAME_1", "NAME_2", "NAME_3")]
result$ID <- ID
print(result)

QUESTION

Now the stuff become tricky because I need to associate additional and deeper information like road_name and road_type.
This information are contained in the shapefiles created on OpenStreetMap and available at: download.geofabrik.de/europe/italy.html.
I loaded the shapefile in R obtaining a SpatialLinesDataFrame:

require(rgdal)
shapefile_roads <- readOGR(dsn = "./road", layer = "roads")

Then, I naively tried to apply the same technique as for joining SpatialPoints and SpatialPolygonsDataFrame:

result <- over(my_spdf, shapefile_roads)

Clearly, the result is just NA. One possible reason that came into my mind was that the coordinates of my_df are not in the exact position of the Lines in shapefile_roads, therefore, I should need some kind of radius parameter. However, I am not really sure.

Can you suggest me the correct approach to perform this spatial join between my SpatialPoints and the attributes of the SpatialLinesDataFrame obtained from the road_shapefile of OpenStreetMap?

Please if something is not very clear do not hesitate to ask.


texts: Split line with UNIQUE NOT NULL attribute

texts:
I have the same problem described here
I am working with QGIS 2.18.10 and postgis 9.4.

Is there a way to make QGIS able to split the line and insert in the constrained field the expression "nextval('fogn_tronchi_cod_tronco_seq'::regclass)"?

My table is like this:

     gid        | integer                | non null preimpostato nextval('fogn_tronchi_gid_seq'::regclass)
 anno       | numeric(10,0)          |
 codice_via | numeric(10,0)          |
 diametro   | numeric(10,0)          |
 id_n1      | character varying(255) |
 id_n2      | character varying(255) |
 materiale  | character varying(255) |
 nome_via   | character varying(255) |
 note       | character varying(255) |
 pendenza   | numeric                |
 pend_ind   | numeric                |
 proprieta  | character varying(255) |
 qfondo_mon | numeric                |
 qfondo_val | numeric                |
 regime     | character varying(255) |
 stato_cons | character varying(255) |
 tipologia  | character varying(255) |
 sollev_rec | character varying(255) |
 recett_rec | character varying(255) |
 comune     | character varying(255) |
 collaudo   | character varying(255) |
 accesso    | character varying(255) |
 prop_serv  | character varying(255) |
 gestito    | character varying(4)   |
 criticita  | character varying(255) |
 the_geom   | geometry               |
 cod_tronco | integer                | non null preimpostato nextval('fogn_tronchi_cod_tronco_seq'::regclass)
Indici:
    "fogn_tronchi_pkey" PRIMARY KEY, btree (gid)
    "fogn_tronchi_cod_tronco_key" UNIQUE CONSTRAINT, btree (cod_tronco)
    "fogn_tronchi_the_geom_gist" gist (the_geom)
Vincoli di controllo:
    "enforce_dims_the_geom" CHECK (st_ndims(the_geom) = 2)
    "enforce_geotype_the_geom" CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL)
    "enforce_srid_the_geom" CHECK (st_srid(the_geom) = 3003)


texts: PostGIS: ST_Clip() raster in SRID:32617 with a polygon

texts:
I'm trying to ST_Clip() raster in `SRID:32617' with a polygon.
What I did

ST_Transform(polygon, 32617): It fails, couldn't project point() exceeded limit (-14)
ST_Transform(rast, 4326, 'NearestNeighbor'): It runs really slow (more than 700,000,000 ms, then I terminated before it finished) about 1300 rasters (tiles)

My polygon got lon/lat from Google map:

ST_GeomFromText('POLYGON((105.78873 21.38659, 106.49185 21.38659, 105.78873 20.71758,
106.49185 20.71758, 105.78873 21.38659))', 4326)

I'm using postGIS 2.1

UPDATE:
My raster is retrieved from Landsat8. I downloaded the raster from earthexplorer.usgs.gov
The raster is imported into postGIS by raster2pgsql with tile-size 500x500, SRID is kept as original (32617,32617, 32632, 32648, 32649, 32658, 32659, 32660))

UPDATE:
Here's Map projection in metadata text of one scene.

GROUP = PROJECTION_PARAMETERS
MAP_PROJECTION = "UTM"
DATUM = "WGS84"
ELLIPSOID = "WGS84"
UTM_ZONE = 49
GRID_CELL_SIZE_PANCHROMATIC = 15.00
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
ORIENTATION = "NORTH_UP"
RESAMPLING_OPTION = "CUBIC_CONVOLUTION"
END_GROUP = PROJECTION_PARAMETERS

UPDATE
The polygon above is invalid. Here's correction:
ST_GeomFromText('POLYGON((105.78873 21.38659, 106.49185 21.38659, 106.49185 20.71758, 105.78873 20.71758, 105.78873 21.38659))', 4326)


texts: Is it possible to only have one Feature Attribute Window open in QGIS?

texts:
I suppose it's by design that when you click through(identify) your features, for every feature it opens the attribute window (which is good, i have enabled the checkbox to automatically open the window).

My question is whether it's possible to limit those, so that the program only shows 1 feature attribute window? Sort of an "auto-close" checkbox?


texts:
I need to import a csv file into a local table. The code I have written is creating the table and populating it, but it is not separating the fields.

My csv file looks like this:

Serial_Num,SB_Lot_Number,Acct_Status,METER_SIZE,Cust_Name,Last_Name,SERVICE_AD,STREET_NAM
11527815,309,Active,5/8",first,last,1111111111,1373,North 10th St
86822741,310,Active,5/8",first,last,2222222222,1367,North 10th St
11527819,312,Active,5/8",first,last,3333333333,1353,North 10th St
11527817,313,Active,5/8",first,last,4444444444,1345,North 10th St

However the output table looks like this (Had to block out some info):

As you can see, instead of creating individual fields for each comma, it is lumping everything together.

Here is my code:

import csv
import os
import arcpy

newFile = r"C:\Users\jraes\Documents\NastyDrankinProjects\NewTestFile.csv"

out_gdb = r"C:\Data\temp.gdb"

arcpy.TableToTable_conversion(newFile, out_gdb, 'tempTable')

It looks like it should be fairly simple. I also did not see in the documentation where I can set a delimiter.


texts: Save user specified settings on flash drive

texts:
With some help in this forum I already managed to successfully setup Quantum GIS 1.8 on a flash drive (see my post on Running QGis from flash drive).

QGis starts and runs stable on different computers, but I still have problems to port also user speficied settings. To port the settings I use the [--configpath] option as suggested in the QGis-Manual. According to this I changed the qgis.bat batch file this way to start QGis with a qgs-project file:

start "Quantum GIS" /B %OSGEO4W_ROOT%\apps\qgis\bin\qgis.exe --configpath %DRV_LTR%\qgisconfig\ %DRV_LTR%\Daten\TESTPROJEKT.qgs

In general this works fine -  QGis starts with the determined project file on different working stations. And installed Plugins, for example, also appear installed when the flash drive is connected to another working station. But it seems, that the [--configpath] option does not port user speficied settings themselves: Installed Plugins do not appear checked when set to another workstation, removed toolbars on one machine appear on the other, so do by „Customization“ hidden parts of QGis.

It seems, that these UI-settings still are saved locally. When removing all entries in the Windows registry, keys appear from new in HKEY_CURRENT_USER\Software\QuantumGIS after running QGis.
I was hoping that the [--configpath] option could also save user specified settings on the flash drive. The overall goal of my work is to supply non-GIS experienced project partners with an easy to handle GIS-Data-Viewer, that can be used without installation. Especially with the new „Customization“ dialog QGis 1.8 seems to be perfect fort that, if not used Menus and Toolbars could just be turned off or hidden.

I am thankfull for any suggestions!


texts: Using Python in the field calculator - if/else over multiple fields

texts:
I'm trying to work out how to populate a field depending on the lowest score across 6 other fields. I've read several other questions relating to if/else statements but I'm only just starting to learn Python so I'm struggling to put together any code that works.

I'm trying to populate a 'Main Failing' field, where for each object ID in the feature class the worst scored field out of 6 other fields will be noted.

The other fields are: 'Interconnectivity', 'PropCover', 'PatchSize', 'ShapeIndex', 'Naturalness' and 'EdgeNat' and they are all short integer, with high numbers being better than low numbers.

So for example if 'Interconnectivity' had the lowest score across the 5 fields for one objectID, I would like to work out how to populate the 'Main Failing' column for that objectID with text saying 'Interconnectivity'.

Edited:
Here's an example of the code I've tried so far with Branco's help. It's running without errors but is either leaving the MainFailing field blank (if 'return resultList[0]') or is populating it with 0's (if return resultList[1]).

Codeblock:

myFunction( !InterconnectivityScore!, !PropCoverScore!, !PatchSizeScore!, !ShapeIndexScore!, !NatScore!, !EdgeNatScore!)

Pre-logic script code:

def myFunction ( interconnectivityScoreValue, propCoverScoreValue, patchSizeScoreValue, shapeIndexScoreValue, natScoreValue, edgeNatScoreValue):
resultList = ["",0]
if interconnectivityScoreValue <= resultList[1]:
    resultList[0] = "Interconnectivity"
    resultList[1] = InterconnectivityScore
if propCoverScoreValue <= resultList[1]:
    resultList[0] = "ProportionalCover"
    resultList[1] = PropCoverScore
if patchSizeScoreValue <= resultList[1]:
    resultList[0] = "PatchSize"
    resultList[1] = PatchSizeScore
if shapeIndexScoreValue <= resultList[1]:
    resultList[0] = "ShapeIndex"
    resultList[1] = ShapeIndexScore
if natScoreValue <= resultList[1]:
    resultList[0] = "Naturalness"
    resultList[1] = NatScore
if edgeNatScoreValue <= resultList[1]:
    resultList[0] = "EdgeNaturalness"
    resultList[1] = EdgeNatScore
return resultList[1]


texts: How to update Trackball controls at three js

texts: I have a project at three.js created from qgis to threejs export where i need to update the trackball controls on window resize.. I need to do that updating the whole controls calling the function with new input variables. Recreating the controls causes crash and i dont want to delete the controls and create new one cause of the garbage this might cause. I have a similar question https://stackoverflow.com/questions/28577347/how-to-update-createcontrolss-function-input-variables-at-javascript but i think this is more complete question.


texts: Perform Radiometric Correction of Landsat images using ArcGIS Desktop

texts:
I need to run NDVI on a study area that is covered by more than one Landsat image.

The NDVI must also be run over numerous years using different Landsat sensors (eg: Landsat 7 Level 1, Landsat 4-5 Level 1 etc).  My research has indicated that normalization across images (temporally, spatially and by sensor) are needed in order to get accurate NDVI values over time.

I believe that radiometric correction will accomplish this however, since I have no experience in remote sensing I do not know how to proceed and I'm looking for step by step instructions.

My understanding is that this should be doable in ArcGIS Desktop.

How do I perform radiometric correction of Landsat images using ArcGIS 10.2 and Spatial Analyst?


texts: Changing scale bar units to feet in PyQGIS?

texts:
I'm trying to change the default QGIS scale bar from meters to feet. From the API documentation it looks like I should be able to use the QgsComposerScaleBar.applyDefaultSize() function.

I thought this would work

scalebar = QgsComposerScaleBar(c)
scalebar.setComposerMap(composerMap)
feet = QgsComposerScaleBar.ScaleBarUnits.Feet
scalebar.applyDefaultSize(feet)
c.addItem(scalebar)

but I get the error

Object::disconnect: No such signal QObject::extentChanged()

Any ideas on how to do this, or what this error means?

I'm using QGIS 2.6.1-Brighton and the CRS for my dataset is EPSG:26917 - NAD83 / UTM zone 17N.


texts: ArcGIS Pro 2.0 and Anaconda environments

texts:
I have had to reinstall ArcGIS Pro 2.0 and have attempted to recreate an Anaconda environment, but am running into some problems. When I create a new environment from the initial ArcGIS Pro environment, I can install packages using pip, but then when I try to load them from within the cloned environment, the modules can't be found. Here are the steps I performed.

Open Command Prompt as provided by ArcPro (not in the GUI)
conda create --clone old --name new
activate new
pip install package
python -c "import package"

This results in an error message. Within the new environment, sys.executable returns the old environment. If I provide the hard path the python executable in the new environment, I can load the package. But this behavior doesn't seem right. I am assuming that when I try to open python within an Anaconda environment it should use the environment's python version.


texts: NULL to 0 values in vactor attribute table in QGIS

texts:
I need to converte NULL values to 0 values with the overall goal of sum up two columns in attribute table in QGIS. Since one of these colums has NULL values, the new column with the sum also appears with NULL value (instead of preserving numeric value of other column).

I have tried Selecting column by all NULL value and putting 0 in Expression Field but that does not work.

Any ideas on how to do it?


texts: Select all polygons (i.e. properties) that can fit a given rectangle

texts:
How can i identify all polygons that are large enough to contain a given rectangle? (Tabular area is not sufficient as i want to exclude polygons that are long and skinny.)

What i want to do is select all properties that are large enough to fit a given #of parking stalls. Some lots are irregularly shaped like flag lots or very narrow. Given the range of shape configurations, i want to be able to select only those lots or polygons that can fit say 4 parking stalls, configured as a minimum bounding rectangle.


texts: Basemap misalignment between Digitalglobe Maps API to Maxar EarthWatch basemap products

texts:
The resulting of features created for Project2018 using DigitalGlobe Maps API (now decomissioned) and Project2019 using Maxar EarthWatch basemap (replacement product) are not aligned.

Is anyone aware of coordinate system parameters of these products and if there is a difference?

There is a 5-10m difference between features as a result of using the 2 basemaps products.


texts: What to use instead of 'Table' input in QGIS 3 modeller

texts:
In QGIS 2 graphical modeler I had an input parameter 'Table' for spreadsheets, csv etc. with no geometry. By now, I'm re-creating all my existing QGIS 2 models from scratch (cp. Using QGIS 2.x processing model in QGIS 3.0?) and miss this input parameter 'Table' in QGIS 3.

What should I use instead?


texts: how is the command of ogr2ogr for converting the adf to shp?

texts:
What is the correct format of ogr2ogr for converting the ADF grid to SHP file? I tried this, but it fails.

$ ogr2ogr -where "OGR_GEOMETRY='Point'" -f "ESRI Shapefile" SHP.shp ADF.adf


texts: Download certain POIs for a whole country and transfer to PostgreSQL database

texts:
I'd like to extract certain data from OSM for a whole country (Germany). The data I'd like to obtain are POI's, but obviously not all of them. I'd like to filter only the supermarkets for example and transfer the directly into a PostgreSQL database.

So far I understood that I can archive that for example with the Overpass API (via Overpass Turbo). I already tried two options:

Overpass-Turbo (browser-based)

QuickOSM-Plugin for QGIS

In Overpass Turbo I filtered through the query-wizard the POIs (shop=supermarket) and set a bounding box fitting Germany inside. It worked only once though (when I chose only nodes and due to access limitations...), but I got also POIS around Germany (as the BBOX doesn't match with the country border). Therefore I tried to use the area-tags (addr:country=DE) to filter only supermarkets in Germany, but then some Points where left out because missing country-tags... anyways when choosing nodes, relations and ways, the query collapses after timeout.

Then I tried the QuickOSM-Plugin for QGIS using the area-code (3600051477) to filter only the supermarkets inside Germany. Same thing – the area seems too big to rum the query successfully. After 180 seconds (max. timeout I guess…) the query aborts. Therefore I tried the same query with sub-regions (like German provinces, for e.g. Berlin, Bayern, etc.) and it worked.

I mean do several queries for each province in a row to fetch the data and merge it after all together, but there should be a more convenient way.

I downloaded the germany osm.pbf file from Geofrabrik and transferred it with osm2pgsql into my PostgreSQL database. It's a huge file (3 GB- takes a lot of disk-space!) and needed several hours to transfer (even with cached 16 GB ram…) and I have now loads of data in my database, which I actually don't need. Altogether not the most satisfying solution as well.

I heard that I could prefilter the .osm data with osmosis before transferring it with osm2pgsql.

Do I have to download the whole dataset .osm.pbf, filter it with osmosis and transfer it with osm2pgsql or is there any other solution, maybe even regarding overpass that I do not know about?

It would be so much more efficient to download only the data I actually need!


texts: Finding featureName in given XY mappoint in ArcMap using ArcObjects?

texts:
I have new in ArcGis Dev. I am trying to create a tool. In the tool when cliked on the feature it will tell you about the XY point of the feature and the name of the feature. I am getting the XY. When I am trying to get the feature Name like if I click on Structure it will show the message "NE.Structure".

Apart structure, all the features I am getting null exception("System.NUllReferenceException occured in the dll but not handled by user code") and 2ndly what I tried is selecting the feature which I think is changing the map state. Can anyone please help. Here is what I tried

public override void OnMouseDown(int Button, int Shift, int X, int Y)
        {
            // TODO:  Add MapControl.OnMouseDown implementation
            screenPoint = new ESRI.ArcGIS.Geometry.Point();
            screenPoint.X = X;
            screenPoint.Y = Y;

            IScreenDisplay screenDisplay = (m_application.Document as IMxDocument).ActiveView.ScreenDisplay;

            mapPoint = screenDisplay.DisplayTransformation.ToMapPoint(X,Y);

            IMxDocument pMxDoc = (IMxDocument)m_application.Document;
            IMap pMap = (IMap)pMxDoc.FocusMap;
            mapPoint.SpatialReference = pMap.SpatialReference;

            IEnvelope env = new EnvelopeClass() as IEnvelope;
            // env.PutCoords(mapPoint.X, mapPoint.Y, mapPoint.X, mapPoint.Y);

            ITopologicalOperator pTopOp = mapPoint as ITopologicalOperator;

            ISpatialFilter pSpatFlt = new SpatialFilterClass();
            pSpatFlt.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
            pSpatFlt.Geometry = pTopOp.Buffer(10);

            ESRI.ArcGIS.esriSystem.UID pFeatLayerUID = new ESRI.ArcGIS.esriSystem.UIDClass();
            pFeatLayerUID.Value = "{40A9E885-5533-11d0-98BE-00805F7CED21}"; /* refer http:// help.arcgis.com/EN/sdk/10.0/ArcObjects_NET/componenthelp/index.html#/Layers_Property/001200000m9t000000/ */
            IEnumLayer pEnumLayer = pMap.get_Layers(pFeatLayerUID, true);
            ILayer pThisLayer = pEnumLayer.Next();

            //IGeometry selectionGeometry = (mapPoint as ITopologicalOperator).Buffer(100); // make a buffer
            pMap.SelectByShape(mapPoint, (m_application as IMxApplication).SelectionEnvironment, false);

            IFeatureLayer pThisFeatLayer = pThisLayer as IFeatureLayer;

            IEnumFeature enumFeatures = (IEnumFeature)pMap.FeatureSelection;

            IFeature thisFeature = enumFeatures.Next();

            do
            {

                IGeometry thisGeom = thisFeature.ShapeCopy;

                   IFeatureCursor pFtCur = pThisFeatLayer.Search(pSpatFlt,true);
                    IFeature pFt = pFtCur.NextFeature();

                    IFeatureClass fC = (IFeatureClass)pFt.Class;

                    MessageBox.Show(fC.AliasName.ToString());

                thisFeature = enumFeatures.Next();
            } while (thisFeature != null);

            MessageBox.Show("X position is " + mapPoint.X.ToString() + "Y position is :" + mapPoint.Y.ToString());
        }
}

The exception I am getting is like :

This above error is coming in this line # IFeatureClass fC = (IFeatureClass)pFt.Class;

You asked me what I mean by changing the map state. In case user already selected any feature or an area, in that if my tool is trying to select the clicked feature, I am afraid it would deselect the user selection.


texts:
I'm new with QGIS and Data Plotly and I have a problem with pie charts feature for Plotly.
I load a shape in QGIS and I use it with Atlas to display a pie chart depending of the selected area on the shape.
In my shape I have a table with cities, Zip Code and a percentage.
I want to display the pie chart with the percentage for the selected city but I don't understand how to do this.

Edit for csk: I downloaded Data Plotly and I need to make a data vis with QGIS. I make a layer for each view and I add a general map and a zoom map in my project. Now I need to display my pie chart for each area I selected. For each area I have in my shape file a colum "CO2 percent" and I need to display this data for each area in my pie chart.
But there are 2 problems: Data Plotly show the percent of the value of my CO2 column on all the values of this column. And not just for example 23 on 100%.
2nd problem is Data Plotly show me in the pie chart all the area around my selected area

I used this tutorial: https://anitagraser.com/author/oneil974/
Some pictures:
This is my table

The rules for the shape

And results without the pie chart with Data Plotly and after with the pie chart

So as you can see in the Table screen I want to display in my pie chart the CO2 value as a percentage. Example for (97360|Apatou|5) I want 5% on 100% on my pie chart and I want to see only 1 "Nom_com" not the "Nom_com" around.


texts:
How can I move the toolbar if its holder is out of screen , like this?

https://i.stack.imgur.com/q6xv9.png


texts: Zoom-in to by coordinates of location in OpenLayers 3

texts:
I'm trying to zoom in to location by coordinates after click on one item from autocomplete list.

    $("#template-custom").easyAutocomplete({
        url: "Toponimos-red.json",
        getValue: function(element) {
            return element.Nombre;
        },
        template: {
            type: "custom",
            method: function(value, item) {
                return  value + " , " + item.Tipo;
            }
        },

        list: {
            match: {
                enabled: true
            },
            maxNumberOfElements: 20,

            onChooseEvent: function() {
                    var selectedItemCoordX = $("#template-custom").getSelectedItemData().CoordX;
                    var selectedItemCoordY = $("#template-custom").getSelectedItemData().CoordY;
                    var valcoordx=parseFloat(selectedItemCoordX);
                    var valcoordy=parseFloat(selectedItemCoordY);

                    var coordinate= [valcoordx, valcoordy];

                    console.log(valcoordx, "-", valcoordy, "-", valcoordx+valcoordy, "-", coordinate);

                    map.getView().setCenter(ol.proj.transform(coordinate, 'EPSG:25830', 'EPSG:3857'));
                    map.getView().setZoom(10);

            },
        },
        theme: "round"
    });

I know that the problem is:

 map.getView().setCenter(ol.proj.transform(coordinate, 'EPSG:25830', 'EPSG:3857'));

Because if I remove this line, the zoom runs.

Why doesn´t it work?  What is the problem?

I have been doing some tests and the problem is the transformation from EPSG:25830 TO EPSG:3857.  If I do the same but the coordinates are in EPSG:3857 it works.

How can I do it with EPSG:25830?

Must I add the plugin proj4?


texts: Creating heat map or density map from stacked polygons in ArcMap?

texts:
I have a polygon shape file that has many overlapping features. Each feature represents an area used by polar bears. So some of the polygons overlap and some don't, and there are separate clusters of overlapping areas.

I would like to get a measure of how much overlap there is in a particular area. But I have no clue how to do this. I'm pretty new to GIS and am still learning the ropes, it just seems I can't find an answer.


texts: How to generate a DEM from a 3D point cloud in MATLAB?

texts:
I'm trying to generate a DEM (or a DSM if that works for this purpose) from a .ply point cloud in MATLAB so that I can generate a hillshade from it and this is what I have managed so far but have no clue where to go from here:
% Input cloud point
[FileName,PathName] = uigetfile({'*.ply';'*.txt';'*.las';'*.*'},'Select orginal point cloud');
ptCloud=pcread([PathName,FileName]);

% Initialize required parameters
points=double(ptCloud.Location);
x=points(:,1);
y=points(:,2);
z=points(:,3);
stepSize = 0.25;

ptMin = min(points);
ptMax = max(points);

minX = ptMin(1);
minY = ptMin(2);
maxX = ptMax(1);
maxY = ptMax(2);

[xq,yq] = meshgrid(minX:stepSize:maxX, minY:stepSize:maxY);
F = scatteredInterpolant(x,y,z,'nearest');
Z=F(xq,yq);

How am I supposed to calculate scale and shift values for this particular point cloud? (I don't know the coordinates system or any other metadata regarding the point cloud).


texts: How to find data orientation atfer an ICA

texts:
I'm trying to cluster tractor sensor data by row in a python script. For the moment I'm trying to find data orientation after a FastICA (sklearn), vertical or horizontal. Because for same data, FastICA can orientate results differently.

Until now, I was searching for an empirical rule with the 2 matrix FastICA().components_ and FastICA().mixing_. But I found nothing that works.

Is somebody have an better idea?

Here is data and my code below:

import geopandas as gpd
import numpy as np

shp_name = 'rawdata379'

shp_source = ('/home/jovyan/scripts/physiocap/Test_Data/{}.shp'.format(shp_name))
source = gpd.read_file(shp_source)
source = source.to_crs({'init' :'epsg:3857'})

src_coord = source[['LONGITUDE', 'LATITUDE']]

from sklearn.decomposition import FastICA
import pandas as pd
from shapely.geometry import Point

rng = np.random.RandomState(42)

pca = FastICA(n_components=2, algorithm='parallel', whiten=True, max_iter=100)
src_coord_pca = pca.fit_transform(src_coord)

src_coord_pca_df = pd.DataFrame(src_coord_pca).rename(columns={0:'X', 1:'Y'})
src_coord_pca_df['geometry'] = list(zip(src_coord_pca_df['X'], src_coord_pca_df['Y']))
src_coord_pca_df['geometry'] = src_coord_pca_df['geometry'].apply(Point)
src_coord_pca_gdf = gpd.GeoDataFrame(src_coord_pca_df, geometry='geometry')
src_coord_pca_gdf.plot(figsize=(9, 9))

compo, feat = pca.components_
print('compo :', compo)
print('feat :', feat)

feat, compo = pca.mixing_
print('compo :', compo)
print('feat :', feat)

#trial of empirical rule
if (pca.components_[0][0] > 0 != pca.components_[0][1] > 0):
    print('vertical')
    src_coord_pca_df = pd.DataFrame(src_coord_pca)[0]

else:
    print('horizontal')
    src_coord_pca_df = pd.DataFrame(src_coord_pca)[1]


texts:
I have GeoJSON multipoint layer and would like to create a line network (GeoJSON multiline) out of it. Best case the points would be snapped to a road. I have only coordinates and no additional information for these points.

I've tried to use QGIS "Points to path" tool, but i have no order for the points, so that's not really the result I'm looking for.
I think the steps would be something like this:

Remove points that are close to each other
Extract street layer from OSM
Snap points to streets
Connect closest points

I don't really know where to start and what tools to use for this.


texts: Enable "Render layers in parallel using many CPU cores" from python console QGIS

texts: How can I enable "Render layers in parallel using many CPU cores" option from QGIS python console?


texts: How can i uncheck leader lines show in piechart using ArcObject?

texts:
I'd like to uncheck the leader lines of my piecharts programmatically . I surf in the internet and found nothing about this.It's simple to uncheck the leader lines in ArcMap. Just go to properties of the piechart and Uncheck Show Leader lines but How can i uncheck the capability using ArcObject? What Interface and Property should i use to unchek th leader lines ?
I'm using ArcGIS 10, ArcEngine 10,c# language
Thanks for Advance.


texts: Append count of unique values to dictionary keys and then sum the occurrences

texts:
I am trying to compare the occurrences of codes from a column within feature classes against a list from a spreadsheet column and then add these to a dictionary where the feature class name is the key and the values are a count of unique occurrences between the two columns. I then need to sum all the value counts for each key. I can make this work to a degree, the problem is that codes are being appended and counted multiple times. How can I make this code work so that code occurrences are only counted once? Here is what I have so far:

FacID = []
with arcpy.da.SearchCursor(IDL,['FACILITY_ID']) as cursor:
    for row in cursor:
        FacID.append(row[0])
FC_FacID = {}
for element in CCMList: ## CCMList is just a list of feature classes (Building, Bridge, etc)
    with arcpy.da.SearchCursor(element,['facilityIDFK']) as cursor_ig:
        for row in cursor_ig:
            key = (element)
            code = (row[0])
            if key not in FC_FacID.keys():
                FC_FacID[key] = []
            if code not in FC_FacID.values():
                FC_FacID[key].append(FacID.count(code))
FC_FacID_Count = {}
for name,code in FC_FacID.items():
    FC_FacID_Count[name] = sum(code)

Sample data:
FacID = [u'NFA001, u'NFA002, u'NFA003]
Feature Class 'facilityIDFK' column: Building(NFA001, NFA001, NFA002, NFA005) #In this example NFA001 appears twice but needs to only be counted once
FC_FacID should look like {u'Building: [1, 1, 0]}
FC_FacID_Count should look like {u'Building: 2}


texts: How to save and load vector layer attribute value map, editable state in .qml style file?

texts:
I just upgraded to Qgis 2.4 and there were problems with style file:

When I load style file for Qgis <= 2.2, the color of features is
right but when I edit them, the value map (value - description pair)
is lost.
The value map I created for layers in Qgis 2.4 is not
kept in .qml style file when I choose 'save style' as 'Qgis layer
style file'.

Is this a bug in Qgis 2.4?

ADDED INFORMATION:

Sample files: http://www.mediafire.com/download/bbk02sj3jl9iwtc/test_data.7z

For line layer line.shp with 3 attributes: id, name, type - I have style file 2_2_line_style.qml generated from Qgis 2.2 which have following lines:

  <edittypes>
    <edittype labelontop="0" editable="0" type="0" name="id"/>
    <edittype labelontop="0" editable="1" type="0" name="name"/>
    <edittype labelontop="0" editable="1" type="3" name="type">
      <valuepair key="first type" value="1"/>
      <valuepair key="second type" value="2"/>
      <valuepair key="third type" value="3"/>
    </edittype>
  </edittypes>

In Qgis 2.2, when I load that line layer with that style and add more lines, there will be a pop up window that allows me to enter values with a drop down list:

But in Qgis 2.4, when I load the same layer with the style file and add (create new) features, the drop down list is gone:

More over, in Qgis 2.4, when I set the value map for that line layer manually, I can get the drop down list but when I save the style as .qml and inspect the file, there isn't any key - value line like the 2.2 style file.


texts: Split Polygon to match length of lines with QGIS

texts:
I have a large polygon shapefile with many features (green), representing waterbody outlines, and a line shapefile (brown) representing river sections. I now want to split the polygon features into smaller parts that match those river sections, but don't know how to do it properly.

I don't mind the software, might be QGIS or ArcGIS, even a Python-script, but let's focus on QGIS for now.  Any suggestions?


texts: Extract longitude and latitude from object in shapefile

texts:
I'm new to GIS and (ESRI) shapefile. I would like to extract GPS-Coordinates or longitude and latitude from every object in a shapefile, mainly from polygons.

Here is a picture with one of these objects:

My plan is to create a GPS-System App where the user can search for that object and my app leading them to the desired object. Is this even possible?

There are somewhere GPS-Coordinates in the Shapefile, because Google Maps can find my shapefile on the map, see here:

So now I would like to know GPS-Coordinates from every Object in that Shapefile, so every colorful object you can see in that picture above.


texts: Can't open Dgn file with CadWorkspaceFactory?

texts:
I am trying to open Dgn File to get data. But I got the error: "Error HRESULT E_FAIL has been returned from a call to a COM component ?" when using OpenFromFile statement, please see my codes:

 Type factoryType = Type.GetTypeFromProgID(
             "esriDataSourcesFile.CadWorkspaceFactory");
 System.Object obj = Activator.CreateInstance(factoryType);
 IWorkspaceFactory factory = obj as IWorkspaceFactory;
 IFeatureWorkspace workspace = (IFeatureWorkspace )factory.OpenFromFile("C:\\DC7.dgn", 0); // Throw Exception

May you help me ?

Regards,
Tai


texts:
I hope to be able to display the coordinates on the map when the mouse moves, not to use marker.
But my code can't run successfully.

This is my code:

<link rel="stylesheet" href="https://leaflet.github.io/Leaflet.markercluster/dist/MarkerCluster.css" />
<link rel="stylesheet" href="https://leaflet.github.io/Leaflet.markercluster/dist/MarkerCluster.Default.css" />


texts: Bulk importing delimited text layers in QGIS Python fail

texts:
I've got a bunch of CSV files to add as text delimited layers to my GIS project. They all have identical formats in the UK EPSG:27700 CRS with a column marked "eastings"(column 0) and another marked "northings" (column 1).

The csv looks a bit like this:

Eastings,Northings,GridRef,Climate,AT,ATA,ATDefault,CT,DAMS,DAMSA,DAMSAdj,MD,SMR,SNR,SiteClass,SMRA,SNRA,SMRN,SNRN,SMRAN,SNRAN,SMRD,SNRD,SNRB,SNRF,SNRG,SMRAspect,Brash,Drainage,Fertiliser,BrashApplied,DrainageApplied,FertiliserApplied,DamsModified,Filter,species,commonname,latinname,modelVersion,modelConfidence,threat,yc,suit,tsuit,limitingFactor,atR,ctR,damsR,mdR,smrR,snrR
363537,473249,SD635732,Medium-High 2050 (A1b/3q0) AWC method,2351,0,2351,7,11,0,11,166,Moderately dry,Poor,Warm - Sheltered - Moist,Moderately dry,Poor,7,2,7,2,0,0,0,0,0,0,No brash present,No drainage installed,No fertiliser,0,0,0,0,All species,AH,Ash,Fraxinus excelsior,3.0,A,Chalara,0,0.0,0.0,SMR,0.99,1.0,1.0,0.9,0.0,0.0
363537,472249,SD635722,Medium-High 2050 (A1b/3q0) AWC method,2366,0,2366,7,9,0,9,168,Moderately dry,Poor,Warm - Sheltered - Moist,Moderately dry,Poor,7,2,7,2,0,0,0,0,0,0,No brash present,No drainage installed,No fertiliser,0,0,0,0,All species,AH,Ash,Fraxinus excelsior,3.0,A,Chalara,0,0.0,0.0,SMR,0.98,1.0,1.0,0.89,0.0,0.0
363537,471249,SD635712,Medium-High 2050 (A1b/3q0) AWC method,2316,0,2316,7,11,0,11,163,Very moist,Medium,Warm - Sheltered - Moist,Very moist,Medium,3,3,3,3,0,0,0,0,0,0,No brash present,No drainage installed,No fertiliser,0,0,0,0,All species,AH,Ash,Fraxinus excelsior,3.0,A,Chalara,8,0.65,0.65,SNR,1.0,1.0,1.0,0.93,0.76,0.65
363537,470249,SD635702,Medium-High 2050 (A1b/3q0) AWC method,2291,0,2291,7,12,0,12,160,Very moist,Medium,Warm - Sheltered - Moist,Very moist,Medium,3,3,3,3,0,0,0,0,0,0,No brash present,No drainage installed,No fertiliser,0,0,0,0,All species,AH,Ash,Fraxinus excelsior,3.0,A,Chalara,8,0.65,0.65,SNR,1.0,1.0,0.99,0.95,0.76,0.65
363537,469249,SD635692,Medium-High 2050 (A1b/3q0) AWC method,2326,0,2326,7,12,0,12,164,Very moist,Medium,Warm - Sheltered - Moist,Very moist,Medium,3,3,3,3,0,0,0,0,0,0,No brash present,No drainage installed,No fertiliser,0,0,0,0,All species,AH,Ash,Fraxinus excelsior,3.0,A,Chalara,8,0.65,0.65,SNR,1.0,1.0,1.0,0.92,0.76,0.65

I have tried importing them using the code listed here, using the editor box of the Python window:

import os.path, glob
from qgis.core import QgsProject, QgsVectorLayer
layers=[]

for file in glob.glob('C:/pathname/*.csv'):
      uri = "file:///" + file + "?delimiter=%s&xField=%s&yField=%s&useHeader=no&crs=epsg:27700" % (",","Eastings","Northings") #(Delimeter, XField, YField) #1st column in file is column 1
      vlayer = QgsVectorLayer(uri, os.path.basename(file), "delimitedtext")
      vlayer.setFieldAlias(0,'X') #these order of columns go: 0, 1, 2 etc
      vlayer.setFieldAlias(1,'Y')
      layers.append(vlayer)

QgsProject.instance().addMapLayers(layers)

but... Literally nothing happens when I run the code.

If I try adding the code to a new script in the toolbox, I get an error Seems there is no valid script in the file.

Can anybody tell me what's going on?


texts: Cannot install GeoPandas on cloned arcgispro env in Anaconda

texts: I have cloned my ArcGIS Pro env into my Anaconda 3 env. My goal was to have access to both the arcpy mod and GeoPandas mod. I tested it in a notebook to make sure arcpy is working and it is. But GeoPandas is not in >conda list. I tried to install it both from the Anaconda Navigator package manager and using >conda install geopandas. Both returned the same error: InvalidSpecError: Invalid spec: =2.7. I've tried to use pip as well. In addition, I tried other solutions, including this solution and this solution where I completely removed and tried to reinstall my gdal package and reinstall geopandas first. It failed and ruined the gdal package in the process. I have tried to clone my arcpgispro-py3 env to see if I can add it that way, but the clone attempt keeps failing. Is there any way to do this or is there some ESRI configuration that makes this impossible?


texts: Change polygon "handedness" for SQL 2008 (reverse polygon vertex order)

texts:
I have a few hundred shapes (polygons and multipolygons) each consisting of tens of thousands of points that I am trying to get into SQL 2008.

Unfortunately, the shapes I have tried to import are "right-handed" (the perimeter of each one is drawn clockwise around the points it contains). SQL server assumes "left-handed" shapes (anti-clockwise around the interior), at least for geography types. This means that SQL assumes that I'm trying to select the entire earth except for my shape. Some people describe this as "inside-out" shapes.

From MSDN, which frustratingly doesn't say which ring orientation one is to use:

  If we use the geography data type to store the spatial instance, we must specify the orientation of the ring and accurately describe the location of the instance.

If you use the wrong ring orientation in SQL 2008, it crashes with the following error (emphasis mine):

  A .NET Framework error occurred during execution of user-defined routine or aggregate "geography":
  Microsoft.SqlServer.Types.GLArgumentException: 24205: The specified input does not represent a valid geography instance because it exceeds a single hemisphere. Each geography instance must fit inside a single hemisphere. A common reason for this error is that a polygon has the wrong ring orientation.

Importing the shapes as geometry instead of geography works fine, but I'd like to use geography if I can.

In SQL 2012, it appears fairly trivial to fix this problem, but I am tied to 2008.

How should I convert the shapes?


texts: Finding total distance of path along points in postgis

texts:
We have a table which contains a list of points, originating from a GPS device's lat/long output every 20 seconds. We want to calculate the total distance travelled within a given time period. We tried out several combination of queries, but it seems that the distance calculated doesn't come near the actual distance travelled.

One of the query that tried so far:

SELECT ST_LengthSpheroid(ST_MakeLine(a.geom), 'SPHEROID["WGS 84",6378137,298.257223563]') AS track_len
    FROM (SELECT geom FROM
gps_data_99 where log_date between '2017-01-21 00:00:00' and '2018-01-22 23:00:00') as a;

But somehow we can't make it work correctly. The total distance travelled today would be around 25kms, and the value returned is:

    track_len
------------------
 203001.154494023
(1 row)

Already looked at:
http://postgis-users.postgis.refractions.narkive.com/JvjGfK5X/newbie-get-total-distance-travelled-from-query-of-points

Total great circle distance along a path of points

We are just starting on PostGIS, so not sure how to get it done correctly.

Table Structure:
CREATE TABLE gps_data_99
(
    gps_data_id bigint NOT NULL DEFAULT nextval('gps_data_99_gps_data_id_seq'::regclass),
    log_date timestamp without time zone DEFAULT now(),
    device_id bigint,
    geom geometry(Geometry,4326),
    CONSTRAINT gps_data_99_pkey PRIMARY KEY (gps_data_id),
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

Environment:
Postgresql 10.1
Postgis: 2.4 (2.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1)
OS: Debian 9.3


texts: Improving accuracy of georeferenced LiDAR data?

texts:
Often I obtain LiDAR data for concept studies and conventional survey (in .dxf format) over a small portion of this. The LiDAR typically covers the broader area and is used for broader studies, catchments, optioneering and concept design. The survey may either be existing survey (for other projects) and provided with the LiDAR, or obtained subsequently specifically for detailed design.

I find that whilst the LiDAR data looks relatively correct, that when compared to the more accurate conventional survey there are often discrepancies. Often the supplier of the data is not familiar with datums or the history of the data, which may also be the cause of some of the discrepancies.

If the discrepancy is relatively uniform, it is relatively straightforward to adjust the LiDAR data. However often this is not the case.

Is there a method to "adjust" LiDAR data to more accurate survey data? I imagine it would be similar to georeferencing  of aerial imagery, but in three dimensions.

I understand that it will never be as accurate, but even if it just looks better I would be happy.


texts: Using PyQGIS to create a new field instead of using the Field Calculator

texts:
I am using QGIS 3.10.
I have imported a csv which contains a column of numeric data (let's call this attribute as x). However, it is imported automatically as strings.
I can easily create a new field by clicking the Field Calculator and typing to_real("x").
Because I need to repeat this many times, I would like to use PyQGIS to do this.
I've been trying various ways but I still can't get it to work. For example, I was following this:
Struggling to use Field Calculator - PyQGIS
Could someone please tell me what I did wrong?
#import data
csv_file = QgsVectorLayer(path_to_csv, 'nameoflayer')
QgsProject.instance().addMapLayer(csv_file)
csv = QgsProject.instance().mapLayersByName('nameoflayer')[0]
#start editing this layer
with edit(csv):
        for feature in selected_feature:
            feature.setAttribute(feature.fieldNameIndex('newfieldname'), to_real(feature['x']))
            updateFeature(feature)

There is no error message but I can't see the new column called 'newfieldname' when I open the attribute table after running the above code.


texts: How to clearly display a huge shapefile with many point features?

texts:
I have two vector shapefiles, each with 54,000 features. When it's zoomed in on a small region, I can clearly see the feature colors:

but as I zoom out, this is what I get:

Any suggestions on how to display the features clearly?


texts:
I have a series of 2D maps with values ranging 0.0 to 1.0 spread all over. The algorithm that produces them is an application of the fractional Brownian motion. Each map measures 101x101 pixels or cells.

The algorithm producing these maps uses a discretized Inverse Fast Fourier Transform (FFT) in the spectral domain. The maps have a varying mean but share a fixed variance which is σ2 = 1. A visual example of such maps is below.

My question: it is said here that applying Moran's I to determine the clustering of a geographic dataset is sensitive to the dataset distribution, especially if the distribution is skewed. Since my algorithm produces maps which distribution is not consistent across all maps, what could be a robust algorithm for determining the clustering of my maps?

The aim is to find a method in which no influence is played by the shape of the distribution on the p-value of the clustering coefficient. The clustering is of course intended to be referred to the red spots in the image below, which are to be seen as my hotspots.


texts:
I'm trying to write a command to see whether a given point lies within a multipolygon boundary (in this case, a watershed). I have a table of watershed boundaries and another table with the point I'm going to use to see which watershed the point lies in.

Would  st_intersects  be the correct command to run? (if so, how would I do that?). If not, what would be the best command to use?

I'm new to postgis/postgresql so any help would be greatly appreciated.

Here are what the tables look like:

Table A has point data:

Table H has polygons:


texts: How to build gdal_calc.py from source in MinGW

texts:
I was wondering how to build the gdal python binding from source using MinGW. In particular, I need to install python utility scripts such as gdal_calc.py.

I've searched around and found this question about Installing GDAL with Python on windows?, but it was about installing gdal binaries. I need to compile from source because I have other libraries in my MinGW/MSYS.

So far, I have downloaded gdal-1.11.2.tar.gz, and I did:

./configure --with-python=/c/Python27/ArcGIS10.2/python.exe
make && make install

/c/Python27/ArcGIS10.2/python.exe above is the python binary from ArcGIS 10.2. I also followed http://trac.osgeo.org/gdal/wiki/BuildingWithMinGW to edit $PYTHON_ROOT\lib\distutils\ccompiler.py.

But I cannot get past the last step make install. After make, I have ./swig/python/scripts/gdal_calc.py in the make folder, but gdal_calc.py is not installed to the msys folder with make install.

I was wondering if anyone can suggest the right way to compile from source here?

The error message I got during make install was:

make[1]: Entering directory `/c/temp/src/gdal-1.11.2/swig'
for dir in python ; do (cd $dir; make install) || exit; done
make[2]: Entering directory `/c/temp/src/gdal-1.11.2/swig/python'
/c/Python27/ArcGIS10.2/python.exe setup.py install --prefix=/usr/local
running install
running build
running build_py
running build_ext
Traceback (most recent call last):
  File "setup.py", line 320, in <module> ext_modules = ext_modules )
  File "c:\Python27\ArcGIS10.2\lib\distutils\core.py", line 152, in setup dist.run_commands()
  File "c:\Python27\ArcGIS10.2\lib\distutils\dist.py", line 953, in run_commands
    self.run_command(cmd)
  ...
  File "setup.py", line 176, in get_gdal_config
    return fetch_config(option)
  File "setup.py", line 129, in fetch_config
    raise gdal_config_error, e""")
  File "<string>", line 4, in <module>
__main__.gdal_config_error: [Error 2] The system cannot find the file specified
make[2]: *** [install] Error 1
make[2]: Leaving directory `/c/temp/src/gdal-1.11.2/swig/python'
make[1]: *** [install] Error 2
make[1]: Leaving directory `/c/temp/src/gdal-1.11.2/swig'
make: *** [install-actions] Error 2


texts: Collision analysis with line diameter in QGIS

texts:
I have some lines representing pipes. They have Z value in the attribute. I need to check if they collide with each other. So I want to know where they intersect/collide in height.
I was wondering if it is possible to do this in QGIS somehow. I now the specified value of the diameter ex 200. Is it possible to make an analysis with the height and the diameter to see where they collide


texts: How can I block WMTS GetCapabilities Link

texts:
I want to block to below link like WMS services.

http://localhost:8080/geoserver/gwc/service/wmts?REQUEST=GetCapabilities

In services, WMS method's has "GetCapabilities" and selected roles only "Admin". When I use this method, User can't display my all layers. I want to same permission for WMTS services.


texts: How to create a polyline in MapInfo Pro by using the snap to nodes function

texts:
I want to create a polyline in MapInfo Pro by using the snap to nodes function, which I can do.

But what I am interested to know is that when a node is added to the polyline, could any of the underlying data attributes/fields of the node that was snapped to (eg ID) be captured into the new Mapinfo table that is being created?


texts: Isolate the name of the city from an address

texts:
I have a data set in CSV format that includes the name of the store, full address, latitude and longitude (from Geocode).

Since address is in full format, is there a way I can get a special column that would only show the city where the store is located?

I don't know how to work in Python/SQL or other program languages.


texts:
I have built a model  in ArcGIS ModelBuilder for which the output is several dbfs having five columns. I want to merge all the output dbfs into a 'master' dbf sheet in the model itself.

How would I go about doing this?


texts: Get Polygon boundaries of all cities or municipalities from whole country in GeoJSON

texts:
Is there a way to get geojson boudaries of all cities in a country ? The only thing I found is to find the geojson of only one city using this:
https://nominatim.openstreetmap.org/search.php?q=Czech+Republic+Prague&polygon_geojson=1&format=json

But I need to generate a geojson file of the whole country(including all small cities etc ..)


texts: Creating duplicate features - many to one - in qgis

texts:
Essentially I have a shapefile of polylines and wish to create duplicates and based on attribute data.

I have route data for individuals as polylines. Then I have attribute data which tells me the number of people and what are the members of the group. The number of individuals is variable ranging  from 1 -10 and I would want to duplicate to create a single line for each individual.

This is to examine footfall. I realise that with a grid I could simply sum the number of members for each polyline togive me footfall. However I wish to examine different members of the group differently and treat their impacts differently, therefore I need single polylines for each member of the group.

A similar question has been asked (below), but this is in ArcGIS and I have only QGIS available.

Creating duplicate features based on many-to-one conversion of related table?


texts: Getting shapefile with change information of two different year's shapefiles?

texts:
I have two shapefiles of different years.

Now, I want to get the change shapefile to show in the form of polygons, and in the table also so that I also can the data of change category


texts: Clipping converts multiband to single band

texts:
I have a multiband raster image. I am trying to extract the AOI by using "Clip raster with polygon" option in SAGA module. The process clips the raster to the exact size of the polygon but the resulting image is in single band and not in multiband as the input. There is no other option in "Clip raster" in SAGA other than input layer and polygon layer. I am using QGIS 3.4.

EDIT:

I tried the option clip raster by mask layer already. I am getting a black rectangle as an output. I used the default setting with match the extent of clipped raster to mask layer extent checked. Both the input and mask layer is having the same CRS.


texts: Converting from Local Coordinate to lat/long and vice verse gives out different results

texts:
I am trying to convert a set of local co ordinates into lat/long and vice versa and I have the code for it taken from this site itself , the issue is When I do these conversions they are not resulting the same.
converting local co ordinates to lat long:
import pyproj
import utm
x = -20509.2462745155
y = -125891.862952247
p = pyproj.Proj("+proj=tmerc +lat_0=34.863361375 +lon_0=138.273186441 +k_0=1.000008632 +a=6378137.0 +b=6356752.314200002 +x_0=-20737.317 +y_0=-126072.022 +units=m +no_defs")
lon, lat = p(x, y, inverse=True)
print(lat,lon)  = 34.864985295652794 138.27568068953587

I also tried the inverse of the above code and it is not working as well.
lat =34.864985295652794
lon= 138.27568068953587
p = pyproj.Proj("+proj=tmerc +lat_0=34.863361375 +lon_0=138.273186441 +k_0=1.000008632 +a=6378137.0 +b=6356752.314200002 +x_0=-20737.317 +y_0=-126072.022 +units=m +no_defs")
nx,ny=p(lat,lon,inverse=True)
print(nx,ny) = 138.50356927081015 36.00068760121446

converting lat long into local coordinates using UTM package.
x, y, zone, ut = utm.from_latlon(34.864985295652794, 138.27568068953587)
print(x,y) = 250960.9300530484 3861456.658975485

as you can see the second piece of code is not exactly returning the X and Y values from the first code snippet.
I am not sure why this is happening , also I can notice that the first couple of values match.
If this is not the proper way to perform this operation then , I would like to know how to approach this.


texts: Checking spatial system reference of shapefile without opening ArcGIS or any other GIS software?

texts: How could you check the spatial reference system of a shapefile without opening up ArcGIS (or any other GIS software)?


texts: How to extract raster values into polygon?

texts:
I have been trying to generate average raster values for a set of overlaying polygons. I would like to put my Raster average value into my overlay polygons. Please help me in this regards.

I am performing these functions on ArcGIS 10.2


texts: Does height difference between beacons results in a more accurate 3D position localisation?

texts:
I'm kinda fresh on the trilateration/triangulation topic.

If you place 4 beacons in each corner of a room at the same height or if you place them on different heights relative to each other will this have an impact on the accuracy of your triangulation/trilateration?

I suppose if you use an error minimizing approach as the Nonlinear least squares it does make a difference as if a closed form algrebraic approach were used, it wouldn't matter?


texts: How to create an Arcade Expression that excludes empty field values

texts:
In ArcGIS Online, I have a point layer with attributes Provider1, Provider2, Provider3.....Provider13. The values are text consisting of people's names.

I would like to configure the popup to concatenate the Providers separated by a semi-colon. However, some points only have a value for Provider1. Other points only have a value for Provider1 and Provider2. Other points only have a value for Provider1, and Provider2, and Provider3... etc.

Currently I simply have:

Concatenate($feature.PROVIDER1, ';', $feature.PROVIDER2, ';',
$feature.PROVIDER3,';', $feature.PROVIDER4,';',$feature.PROVIDER5)

But for points that have only Provider1, the popup shows:

  John J. Doe; ; ; ;

Is there a way to to have the expression "stop" trying to concatenate when it encounters an empty Provider value?


texts: How to convert a QGIS Composer Grid to polygons

texts:
I really like the cartographic aesthetics of the QGIS 3.10 Composer Grid, especially the border created by using the "Zebra" style and "Fit Segment Length" option.

Because of its design, I'd like to use this Composer Grid as the underlying groundwork to generate a points of interest locator table.

Is it possible to export the Composer Grid to a polygon layer so that I can run a point-on-polygon spatial overlay?


texts: How to plot street addresses and streets using Python/GIS?

texts:
I asked this on stackoverflow but did not get a really clear answer.
I am not really familiar with GIS systems in python. I am looking for a library to do the following:

Lets say I have a pandas data frame with one column of street addresses and another of labeled categorical or numerical data. Something like the following for San Francisco:

df=pd.DataFrame({'Location': ['1 Dr Carlton B Goodlett Place','1390 Market Street','1100 Fillmore St' ],
                 'Type': ['G', 'G','C']})

And I want to plot them on a map of San Francisco with the 'Type' column being in red for the data in 'G' and blue for 'C'. How would I a simple plot like this? How would I also include streets on the map? Can an simple example be provided?

I found a library called Geopy to change the data to longitude & latitude if it helps to convert the data. Also read that matplotlib can do this with basemap but it is for areas larger than cities.

Is there a quick way I can do this? Would I have to use ARCGIS?


texts: Nominatim - from Point to Polygon (town, city boundary)

texts:
I'm struggling with Nominatim's Search and Reverse.

What I want, is to reverse geocode a Point and get a (Multi)Polygon. An administrative boundary (town, city scale preferred).

reverse geocoding: sometimes identifies neighbouring town/city
instead.
search: input coordinates. Still identifies neighbouring
town/city instead.

Example:
Nominatim Reverse lookup for Point

Result mentions Bad Berleburg. If you look that up, the point is not in the polygon for Bad Berleburg. Somehow I should derive a boundary for Bad Laasphe I guess (so from the point I actually need to derive this boundary, for which I don't know the name in advance without reverse geocoding).

I'm looking to accomplish this with a free API, preferably Nominatim. This is my fallback method. I'm using a database with GADM shapes for the 2 most common countries, but I don't want to import all the cities in the world in that database by default.

Update: I've been pointed to the Overpass API for this, looking into it now. Looks promising.


texts: Try to load GeoJSON object with OpenLayers 3

texts:
I'm a bit new to working with OpenLayers 3. I'm trying to load a GeoJSON object in OpenLayers 3 but it display the center only point[0,0], which not even in my GeoJSON object, where I'm going wrong.

The code is below:

//CREATE A BASE LAYER
        var raster = new ol.layer.Tile({
          source: new ol.source.OSM()
        });

        // create a source from the feature
        var source = new ol.source.GeoJSON(

            ({
                    object:
                        {
                            'type': 'FeatureCollection',
                            'crs': {
                                'type': 'name',
                                'properties': {'name': 'EPSG:3857'}
                            },
                            'features': [
                                {
                                    'type': 'Feature',
                                    //'properties': {'id': 63},
                                    'geometry': {
                                        'type': 'Point',
                                        'coordinates': [12, 12]
                                    }
                                },
                                {
                                    'type': 'Feature',
                                    'properties': {'id': 62},
                                    'geometry': {
                                        'type': 'Point',
                                        'coordinates': [35.0, 1.0]}
                                },
                                {
                                    'type': 'Feature',
                                    'properties': {'id': 61},
                                    'geometry': {
                                        'type': 'Point',
                                        'coordinates': [34.0, 0.0]}
                                },
                                {
                                    'type': 'Feature',
                                    'properties': {
                                        'id': 56
                                    },
                                    'geometry': {
                                        'type': 'Point',
                                        'coordinates': [33.0, 33.0]
                                    }
                                }
                            ]
                        }
            })
        );

       // CREATE THE LAYER WITH THE DATA
        var vectorLayer = new ol.layer.Vector({
            source: source
        });

       // create an openlayers map object
        var map = new ol.Map({
            target: attrs.id,
            layers: [raster,vectorLayer],
            view: new ol.View({
                //projection:'EPSG:4326',
                center:[0,0],
                zoom:2
            })
        });
    }


texts: Deleting only specific LiDAR points in QTM with Mensuration tool

texts:
I'm using Quick Terrain Modeler (QTM) 64-bit v8.0.3.4 - I have a set of LiDAR data that has some points that I want to delete, but due to their location it is hard (if not impossible) to delete them by merely using the Select Polygon (Screen) tool without deleting other desired points.

So, the way around it is to use the Start Mensuration tool and then select Profile Analysis Tool in order to cut away specific points. However, this can be a bit tedious and not always a one-shot process. I found that while in the Profile Analysis you can select Mask to Area in 3D, in which (temporarily) only the selected area of interest is shown. Once in this view, I choose the Select Polygon (Screen) tool again and select the bad points and then cut them. However, this cut still deletes the points behind the bad ones, even though I've temporarily made everything else disappear.

So my question is this:

Is there a way use the Select Polygon (Screen) and Cut tools in conjunction with the 'Mask to Area in 3D' selection, without cutting/deleting points that are not visible/in the area of interest?

Here's an example of what I'm talking about. The first picture is with Mask to Area in 3D turned on and I'm cutting the faulty points. The second picture is showing the unintended cut, after turning off Mask to Area in 3D.


texts: Does anyone knows If there is a way to calculate Geary's index of spatial autocorrelation on open source gis?


texts: Apply restrictions to publishers on ArcGIS Server?

texts: We have added a specific AD group to the  publisher role on our ArcGIS Server, which is version 10.1.  We also created a folder that has its security settings restricted to this same AD group.  We would like this group to have the ability to publish services, but restrict their ability to publish to just this single folder that their AD group has access to.  Unfortunately, it seems that once a group is added to the publisher role, they are able to access all services on the server and publish in any folder.  Is there a way to restrict where a publisher group can publish services to a single folder and also restrict their access to just the services in that folder?


texts:
I have a PostGis database with a multipolygon geography (as opposed to geometry) column. I'd like to fill it with the country border data available from http://thematicmapping.org/downloads/world_borders.php The query I tried is

insert into "countries"
select "un", "iso3", "name", "region", "subregion",
 ST_GeogFromText( 'POINT(' || "lon"::varchar || ' ' || "lat"::varchar || ')' ),
 ST_Transform( "the_geom", 4326 )
from "world";

Where "world" contains an import of the raw (geometry) data from the link. However, it complains about unknown input geometry SRID -1. I have tried explicitly setting it to 900913 via ST_SetSrid, but that didn't change anything (and I'm not even sure it's the correct SRID to use, just guessing).

So my question is: How do I get a country border data set into my geography column?


texts: Polygon Selection in Manual Editing (eCognition)

texts: eCognition software allows manual editing on segments with the options such as"polygon selection, line selection, rectangle selection and cut object manually". In my study, polygon selection option does not work correctly when determining a boundary and just emphasizes the current segment. So, a new edited closed area cannot be created manually but "cut object manually" works correctly. Is there any setting to overcome this problem?


texts: polygon in polygon in polygon(or more...)

texts:
I have a buildings feature class which contains a building in a building (4 times), and I want to extract each building to a new feature class which  created. I have written a Python code that does half what I want.

it creates a new field (count in)
it makes a feature class
selection by location and then by attribute
creates a new feature class for results
count

I want that script extract the single building and give them the number 1 in "count_ in" field, and then copy the row to the new fc(fc_1).
then extract the building in building and give the number 2 to the feature class which selected, and copy to the new fc(fc_2).
exactly to the other feature (3 bldg and 4 bldg).

this is the code:

#this script selected(by location)bldg in bldg, and print the result by 1 or 2
#or 3 (or more)in a new field("count_in")if was select bldg in bldg.

print "Start"

import time
import arcpy
from time import strftime

arcpy.env.overwriteOutput = 1

FC_BLDG = r"D:\AVI_DAVID\DB_CORRECTION_TEST\test.gdb\poly_test"
gdb_name = "bldg_in_bldg_check"+'_'+strftime("%Y-%m-%d_%H%M")#
GDB_Destination = r"D:\AVI_DAVID\DB_CORRECTION_TEST"#
GDB = arcpy.CreateFileGDB_management(GDB_Destination,gdb_name+'_'+strftime("%Y-%m-%d_%H%M")) #create gdb+location

count = arcpy.GetCount_management(FC_BLDG)
for con in count:
        print con ," TOTAL bldg COUNT"
        print""
        print""

##fieldList = arcpy.ListFields(FC_BLDG)
##fieldListName=[]
##for field in fieldList:
##    fieldListName.append(field.name)
##    if "count_in" not in fieldListName:
##        arcpy.AddField_management(FC_BLDG,"count_in")
##        arcpy.CalculateField_management(FC_BLDG,"count_in","1")

def count_inside(FC_BLDG,c):

    fieldList = arcpy.ListFields(FC_BLDG)
    fieldListName=[]
    for field in fieldList:
        fieldListName.append(field.name)
    if "count_in" not in fieldListName:
        arcpy.AddField_management(FC_BLDG,"count_in")

    if c == 1:
        exp = '"count_in" IS NULL'
    else:
        exp = '"count_in" IS NOT NULL'

    lyr1= arcpy.MakeFeatureLayer_management(FC_BLDG,"lyr1",exp)
    lyr2= arcpy.MakeFeatureLayer_management(FC_BLDG,"lyr2",exp)

    print "count features before selection",arcpy.GetCount_management(lyr1)

    arcpy.SelectLayerByLocation_management(lyr1,"COMPLETELY_WITHIN",lyr2, "","NEW_SELECTION")

    NEW = arcpy.CreateFeatureclass_management(GDB,"BLDG"+str(c),"POLYGON",FC_BLDG)
    arcpy.CopyFeatures_management(lyr1,NEW)

    try:
        arcpy.CalculateField_management(lyr1,"count_in",c)
    except:
        print arcpy.GetMessages(2)

    print "count features after selection",arcpy.GetCount_management(lyr1)
    return arcpy.GetCount_management(lyr1)

c=1
count_inside(FC_BLDG,c)

c=2
count_inside(FC_BLDG,c)

c=3
count_inside(FC_BLDG,c)

c=4
count_inside(FC_BLDG,c)

print "-------OK-------"


texts:
In my android application I add a new feature to a feature layer after processing data in various activities, feature layer is updated but when I return to the main activity to display updated feature layer neither map nor feature layer is displayed In  the log cat I found
com.esri.core.io.EsriServiceException: Request Time-out


texts:
I'm working on a project that routes various modes of transport to one location (a school) in Arc Map 10. I have about 400 routes, so there is a significant amount of overlap as they near the school. Just for simplicity, let's say that 50 students walked, 150 drove and 200 took the bus.

What I'd like to is show the difference in modes of transport by widening the polylines as they overlap. And then, as they get close to 1 or 2 blocks from the school, have them all join. The line for students that walked to school would be 12.5%, students who drove would be 37.5%, and the students who bused would be 50%.

Can you think of a way I would be able to symbolize the data to show this? It's a really specific problem I'm working on, so I've been having trouble coming up with the right solution. What I've done so far is separate the polylines that intersect, but that's where I'm stuck.

Any help would be greatly appreciated!

Thank you!


texts: Most effective way to find spatial relationships in python

texts:
I am trying to find the spatial relationship between two polygon geoDFs as follows:

Find geoDF1.geom within geoDF2.geom else (not within)
Find geoDF1.geom with largest intersection fromgeoDF2.geom else (no intersection)
Find geoDF1.geom with smallest distance from geoDF2.geom
Update geoDF1 with selected column from from geoDF2

I wrote up a function which seems to works fine, but I am just not sure this is the most effective way to do it. I am especially worried about step #3 which right now checks the distance from each one of geoDF2 rows, which can become memory intensive if it is a large table. Perhaps I should create a buffer that limits the distance checked?

import geopandas as gpd

### FUNCTION that checks spatial relationship between two geodata frames - df1, df2:
# Order of operations: Within -> largest intersection -> smallest distance
# Checks relationship of df1.geom (left) to df2.geom
# Updates df1 with new col from df2:: df2_col
# df1_id = unique id

def relate(df1, df1_id, df2, df2_col):

    # join df1_col for all df1_geom within df2_geom
    join_df = gpd.sjoin(df1, df2, how='left', op="within")
    ##sort by col_1
    join_df = join_df.sort_values(df1_id)

    # get rows that are not within
    remaining_df = join_df.loc[join_df[df2_col].isna()].reset_index(drop=True).drop(columns='index_right')

    ##checks to see relations of df1 to df2
    intersects = pd.DataFrame(columns=['id', 'type'])
    # check for each geom1 not within geom2
    for df1_geom, id in zip(remaining_df.geom, remaining_df[df1_id]):
        size = 0
        for df2_geom, type in zip(df2.geom, df2[df2_col]):
            # see size of intersection , if larger (or first) than previous intersection, set 'sel_type' to df2_col to associate with df2 (largest intersection)
            try:
                new_size = df1_geom.intersection(df2_geom).area
                if new_size > size:
                    size = new_size
                    sel_type = type
                    # print("intersects" + str(id) + sug)
            except:
                print("couldn't intersect for df1 with id: " + str(id))
        # if intersection exists add to lists of df1_id, df2_col
        if size > 0:
            intersects = intersects.append({'id': id, 'type': sel_type}, ignore_index=True)

    # sort by df1_id
    ints = intersects.sort_values('id').set_index('id', drop=True)

    # set type for intersecting geoms
    join_df.loc[join_df[df1_id].isin(ints.index), df2_col] = list(ints['type'])

    # #get df1 geom, who are not within and do not intersect with df2_geom
    remaining_df2 = join_df.loc[join_df[df2_col].isna()].reset_index(drop=True).drop(columns='index_right')

    ##!!## JOIN based on distance
    distance = pd.DataFrame(columns=['id', 'type'])
    # check for each df1 with no df2_col
    for geom, id in zip(remaining_df2.geom, remaining_df2[df1_id]):
        size = None
        for df2_geom, type in zip(df2.geom, df2[df2_col]):
            # see size of distance from df2_geom, if smaller (or first) than previous distance, set 'sel_type' to associate with df1_geom (minimal distance)
            try:
                insize = geom.distance(df2_geom)
                if size == None or insize < size:
                    size = insize
                    sel_type = type
            except:
                print("couldn't calculate distance for df1 with id: " + str(id))
        if size != None:
            distance = distance.append({'id': id, 'type': sel_type}, ignore_index=True)

        # sort by df1 id
    dis = distance.sort_values('id').set_index('id', drop=True)

    # set df2_col type for remaining df1 geoms by minimal distance
    join_df.loc[join_df[df1_id].isin(dis.index), df2_col] = list(dis['type'])

    ##UPDATE columns on original df1 with col to join from df2
    df1[df2_col] = [join_df.loc[join_df[df1_id==x], df2_col] for x in df1[df1_id]]


texts: Why won't this simple raster layer in my mapfile show up?

texts:
I am trying to get a simple raster layer to work on MapServer. I am not exactly going to be using geographic data, but I want to "georeference" multiple rasters together to create a map-like interface. Unfortunately, I can't seem to get over this initial step of getting a single raster to show up using my mapfile.

Mapfile:

MAP
    IMAGETYPE   PNG
    EXTENT      -180 -90 180 90
    SIZE        400 300
    IMAGECOLOR  255 255 255
    UNITS       DD
    STATUS      ON
    DEBUG       ON

    LAYER
        NAME    "base"
        DATA    "Base.png"
        STATUS  ON
        TYPE    RASTER
    END
END

My directory structure is like this:

/home/username/public_html/cgi-bin/mapserv.cgi

/home/username/public_html/dir/Base.png

/home/username/public_html/dir/mapfile.map

/home/username/public_html/dir/tmp/

When I visit the map using something like the following link, all i get is a blank image of the correct size but not the actual Base.png image. No pink tiles. No error message. What am I doing wrong?

*http://myurl.com/cgi-bin/mapserv.cgi?map=/home/username/public_html/dir/mapfile.map&layer=base&mode=map


texts: Calculating the distance between Centroid and Point

texts: I am a college student trying to finish a project in ArcMap; I have no professional training and am at a loss.  I need to calculate the distance, in miles, from the centroid of polygons to a specific point.  I have made this point a shapefile.  I think I should be using Python, but I am not sure what expression to use.


texts: Adding legend for specific ranges not for discrete values in QGIS

texts:
I want to add legend for a QGIS proportional symbol map rangewise. That is my data ranges from 57000-630000. So I want the legend as:
57000-60000,
60500-70000,
74000-79000,
81000-87000,
90000-94000,
115000-124000,
125000-226000,
250000-290000,
293000-300000,
350000-373000,
374000-390000,
600000-630000

But in QGIS I am only able to add the legend for a particular value like 57000, 69000. Is there any way to add the legend for ranges?
On adding the legend it displays this:

Layer panel screenshot:


texts: How to split polygons by overlapping feature?

texts:
How can I divide the multiple overlapping polygons into multiple non-overlapping polygons? Here is an example: I would like to divide the long yellow polygon in the image into 3 polygons based on where it touches the other two.

library(spData)
library(sf)

seine_buff_5km = st_buffer(seine, dist = 5000)

tmap_mode("plot")
tm_shape(seine_buff_5km) +
  tm_fill( col = "name", alpha =.5) +
tm_shape(seine) +
  tm_lines()

Of course, I would like to do that for a larger data set with many more polygons. I guess one way to do it is subtract all polygons in an sf object from each other.


texts:
This is a duplicate asked on AutoDesk Forum.

If I form a request like this:

http://usalvwdgis1/mapserver2013/mapagent/mapagent.fcgi?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=2&FID=47403&TYPENAME=Water:WA_LARGE_VALVE

I get a valid response:

<gml:featureMember>
  <Water:WA_LARGE_VALVE>
    <Water:FID>78935</Water:FID>
    <Water:Geometry>
      <gml:Point>
        <gml:pos>498289.041000 5459452.306000</gml:pos>
      </gml:Point>
    </Water:Geometry>
  </Water:WA_LARGE_VALVE>
</gml:featureMember>

But if I add a srsName to the request I get the exact same response:

http://usalvwdgis1/mapserver2013/mapagent/mapagent.fcgi?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=2&TYPENAME=Water:WA_LARGE_VALVE&srsName=EPSG:26910
http://usalvwdgis1/mapserver2013/mapagent/mapagent.fcgi?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=2&TYPENAME=Water:WA_LARGE_VALVE&SRSNAME=EPSG:4326

But I expect to see the results in the requests coordinate system.  Is my request incorrectly formed (maybe add 'urn:x-ogc:def:crs' before the code)?


texts: Is their the effective layer limit for a Geoserver instance running on Ubuntu Linux 16.04?  I assume that each layer translates to a Postgres table, and while I have read that there is no real limit by Postgres on number tables per instance, there are some limits imposed by OS and performance considerations. What count of layers in GeoServer might you begin to experience performance problems?


texts:
I have a table similar to the following attached to geometry destined for inclusion in an esri geodatabase:

What I am trying to do in FME, is create a new field containing a hyperlink to reference the file name and location of the reference number at a fixed networked location. For example:

Is this possible in FME and if so, what would be the best way to go about it? There are occasions where multiple concatenated reference fields exist, so I wonder if it's possible to also concatenate the hyperlinks or if they would be required to be split into separate fields?


texts: Sorting grid cells in fishnet grid by location and number using QGIS?

texts:
I am using QGIS 3.2.0 on a Mac (and also have QGIS 2), and I have a fishnet grid (probably created in Arc) that comprises of nearly 2 million grid cells. As a result, every change I make to the attribute table takes hours, so I can't experiment as I would usually.

I need to label each grid cell such that the north-westernmost cell is "1", and the numbers proceed by rows, until the last number goes into the south-easternmost cell. Using the row-number doesn't work, as that's not how the table is sorted. I've done this before, but cannot remember how.


texts: gdal_translate VRT to TIF has sub-pixel shift

texts:
Problem:
While combining 2 tif images using gdalbuildvrt and gdal_translate, the output tif has sub-pixel shift. (GDAL 3.1.1, released 2020/06/22)
File names:
input files: input1.tif, input2.tif
input files list: input_list.txt
output VRT file: merged_output.vrt
output TIFF file: merged_output.tif

Dropbox links to the two test input files are given below.
Code (two steps):
gdalbuildvrt -srcnodata 0 merged_output.vrt -input_file_list input_list.txt -a_srs EPSG:3413
gdal_translate -of GTiff -r lanczos --config GDAL_VRT_ENABLE_PYTHON YES -ot Byte -co COMPRESS=DEFLATE merged_output.vrt merged_output.tif

Output

Link to test data:
input1.tif input2.tif
Can someone experienced, please share his/her understanding here. I would like to use buildvrt + gdal_translate and not gdal_merge.py due to other concerns.


texts: QGIS 3.10 Attribute form "default value" with current field values, but NOT update automatically

texts:
I want to have a default value from a current field (from another layer, but it is the same, when I use a current field from the current layer).  And I DON'T want to update it automatically! After getting the current value as default value, I want to have the possibility to change the field value individual.

In 3.10 it is no more possible, there applies a warning “to use fields in default value is only working if you use update automatically default value”. But I don't want to update it automativally!:

In 3.4 is was possible and worked fine, I could have a default value from a current field without update it automatically:

Will it work again like in 3.4 in a newer version? Because of the warning it seems to be a wanted change and it doesn't seems to be a bug. But for me it is like a bug, I would like to work with 3.10 or newer. 3.4 I installed for test on an not used computer...


texts:
I have 7 animated gifs (5 animations in each gif) with bounding coordinates of real time radar data. The number of these will grow significantly over the next year or so.

We are currently serving each gif as a separate raster layer using OpenLayers but this is causing problems on slower connections as OpenLayers appears to load the gifs even if not being displayed in the viewport.

I would like to be able to combine and serve all the animations via GeoServer as one layer, probably animated using the GeoServer Animation feature.

This is totally rocket science to me so I am reaching out for suggestions from more experienced users.

I guess we have to break the gifs apart for each animation, combine them into one raster for each time slice (they are all synced), load each georeferenced raster into PostgreSQL then use GeoServer's animation feature to display the layer.

Am open to any and all suggestions.


texts: Calculating the distance between the furthest vertex of a Voronoi polygon and the point from which it was generated

texts:
I need to calculate the distance between the furthest vertex of a Voronoi polygon and the point that generated it. I need to calculate this distance for all Voronoi polygons.

Is there a way to do it automatically using QGIS?

I attach an image to give the idea, in which I would like to calculate the length of the highlighted yellow line. I also report some sample data from my layer containing the Voronoi polygons (TOPOIETI + IMP) and from the layer of the points from which they were generated (NODI TOPOIETI).


texts: QGIS 3.4.3 Refactor fields changes symbols in fields (cp1250)

texts:
I've been using QGIS Refactor fields to change type, re-order and rename a substantial ammount of columns (statistical data). Most of them hold purely numerical values and are as such not material to whatever decoding/encoding happens between various codepages, except for one. This column holds names (place names) which contain characters č, š and ž. Original data I procured is encoded in utf-8. For various practical reasons I exported them as cp1250 (default for all the computers involved) which produced exactly the same file with properly displayed symbols (čšž).
When I applied Refactor fields everything was copied to a new file but č symbol was exchanged for c, while š and ž remained as they were.

I don't see any way to specifiy encoding in Refactor fields and since this is based on Python (which I know defaults to utf-8 in 3.x) means that I have encountered this kind of problem before. In Python that could be corrected by specifiying a codpage (in most cases).

What else could be the source of this issue and, in particular, can it be avoided using the command line version of the Refactor fields tool?


texts: QGIS - Simplification process creates a new polygon - how to avoid it?

texts:
I am trying to simplify a layer made of polygons.
I have used different methods but the new simplified layer has a polygon which intersect with itself, thus creating a "new" polygon.
Please see screenshot below for an example.
The darker line is the new layer and instead of following the original line (you can see all its vertices as red crosses) it cut across to reach the top left vertex.
Is there any method/settings that I can use to force the simplification to follow more strictly the original layer boundary?

I have tried douglas_reduction with a Maximal Tolerance Value of 0.000001 and Percentage of point of 99 and it worked if applied to ONLY one polygon. If applied to the whole layer I still have the result shown in the screenshot.


texts: Downloading and processing raster files in Python?

texts:
I am fairly new to python and seek guidance for a question which might sound trivial to many.

Is there a way to use 'wget' in a python script to download raster files from a server and process them in the same script?


texts:
I have a simple arcpy script to update a field in a point shapefile with info from the polygon feature that it is within.  It takes 9 minutes to do 100 points in arcpy but a spatial join in arcmap is instantaneous. I'm sure there is a quick established way to solve this problem.  Can someone point me in the right direction?

import took 0:00:07.085000
extent took 0:00:05.991000
one pt loop took 0:00:03.780000
one pt loop took 0:00:03.850000
one pt loop took 0:00:03.791000

import datetime
t1 = datetime.datetime.now()
import arcpy
t2 = datetime.datetime.now()
print "import took %s" %  ( t2-t1)
#set up environment
arcpy.env.workspace = "data\\"
arcpy.env.overwriteOutput = True

desc = arcpy.Describe("parcels.shp")
ext = desc.Extent
extent = (ext.XMin,ext.XMax,ext.YMin,ext.YMax)
t3 = datetime.datetime.now()
print "extent took %s" %  (t3 -t2)
fc = arcpy.CreateRandomPoints_management("", "malls.shp", "", ext, 100, "", "POINT", "")
arcpy.AddField_management("malls.shp", 'ParcelID', 'LONG')

rows = arcpy.UpdateCursor('malls.shp',"","",'ParcelID')
for row in rows:
    t4 = datetime.datetime.now()
    pt = row.Shape.getPart()
    for polyrow in arcpy.SearchCursor('parcels.shp'):
        t6 = datetime.datetime.now()
        poly = polyrow.getValue('Shape')
        if extent[0]<pt.X<extent[1] and extent[2]<pt.Y<extent[3]:
            if poly.contains(pt):
                print "works"
                row.ParcelID = polyrow.Parcels_ID
                rows.updateRow(row)
                break #we can stop looking for matches since
        t7 = datetime.datetime.now()
        "a full poly loop took %s" % (t7-t6)
    t5 = datetime.datetime.now()
    print "one pt loop took %s" % (t5-t4)

print datetime.datetime.now() -t1


texts: Unable to import QGIS after setting PYTHONPATH and LD_LIBRARY_PATH on Mac

texts:
I've set my environment variables on Mac to:
export PYTHONPATH="/Applications/Qgis.app/Contents/Resources/python"
export LD_LIBRARY_PATH="/Applications/Qgis.app/Contents/Resources/lib"

but am still getting
Traceback (most recent call last):
    File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'qgis'

Pretty new to Mac (coming from Linux) so wondering if I might be setting the wrong environment variable, but couldn't find any other alternatives on here.


texts: What map library does this census visualization use?

texts:
Link: http://projects.nytimes.com/census/2010/explorer

How can I also render that specific kind of map? Does it come from some popular library?

It doesn't look like Natural Earth or drawn with Mike Bostock's TopoJSONs. It says its from Google at the bottom right, but then why does it look different from maps.google.com?


texts: Calculating distance between coordinates from point to a box?

texts:
I have two different WGS84 coordinates, A and B. A is just a (lat,lon,z) point and B is a (lat,lon,z) point inside the center of a box with the dimension height, width, depth (in meters).

I need to calculate the distance from A to the box, not to the center of the box but to the closest side/surface of the box. Any ideas?

I already have a method to calculate between two coordinates ((lat,lon) to (lat,lon)).

Currently using PostGIS, but open to suggestions for other software.


texts: Definition queries in Python Toolbox tools

texts:
I'm trying to set a definition query within a Python Toolbox tool in ArcGIS Desktop. The code below works when I run it in the ArcGIS Python window, but does not perform the definition query or provide any error messages when I run the script tool that this code is part of. How do I set a definition query up to work in a script tool? This code (along with other code) is within the def execute(self, parameters, messages): block of the toolbox template that ESRI provides (https://desktop.arcgis.com/en/arcmap/10.3/analyze/creating-tools/a-template-for-python-toolboxes.htm)

import arcpy.mapping as mp
        mxd = mp.MapDocument('CURRENT')
        queryValue1 = 2
        queryValue2 = 1000
        lyrs = mp.ListLayers (mxd)
        for l in lyrs:
            if l.name.endswith('rockfall'):
                l.definitionQuery = "{} = {}".format(arcpy.AddFieldDelimiters(lyrs, 'gridcode'), queryValue1) + ' AND ' + "{} > {}".format(arcpy.AddFieldDelimiters(lyrs,'Shape_Area'), queryValue2)
        arcpy.RefreshActiveView()
        arcpy.RefreshTOC()


texts:
I'm a little new to this, is there a way to have features on a layer render only if they overlap with features on another layer?

I have roads of a city, and different areas in the city, and I'd like to only render roads that are entirely within the bounds of the city. I guess that means using a rule-based render on the features in the road layer, I just don't know the syntax of the condition "Overlaps with city layer."


texts: parcels do not display at some scale in SAEx

texts: We have a geodatabase of cadastral parcels of every district. In some of the database, the map is black in some areas, on zooming in/out those areas display whereas other areas become blank. we have tried importing the database to the newer version as well but the problem still exists.


texts: Using multipolygon shape layers in profile tool

texts:
I am working with QGIS for some time. The profile tool is indeed very useful to visualize the 3D problems. My problem now is associated with this tool.
I have a multi polygon shape layer which I assigned z levels by set z option in processing tool hence I believe z levels are assigned to the vertices of the multi polygons I have. I wanted to use this shape layer to plot a cross-section by using the profile tool however got an error as "unprofilable layer" or such. It is understandable since this layer is not a surface or a point data neither a raster. I cannot export this layer as DXF and turn it into a surface since I need the attributes assigned to these multipolygons for my project.
Is there any method so that I can transform this multipolygon shape layer into a kind of surface and/or raster or something else which I can both use it in the profile tool and also keep the attributes assigned?


texts: I want to insert hospital points inside district Polygon. District is divided into 9 polygon called WARD_NO and hospitals have locations according to polygon 1,2,3,4,5..... not according to x and y co-ordinate,so how can  I do so by using .csv file format which have column HOSPITAL_NAME and WARD_NO only. I want to locate those hospitals according to WARD_NO in the respective WARD_NO


texts: Out of bound error updating recently created attribute with memory layer

texts:
I'm creating a new memory vector layer (polygon) copying data from a selected input layer on my plugin (QGIS version 2.18.2 on mac). New data columns (attributes) must be added with results of calculations. (Pseudocode)
    sourceLayer = self.layers[self.confirmedLayerIndex]  #layer used as input
    sourceFeats = [feat for feat in sourceLayer.getFeatures()]
    sourceGeometryType = ['Point','Line','Polygon'][sourceLayer.geometryType()]
    sourceCRS = sourceLayer.crs().authid()

    name = "layername"
    data = matrix with shape (d, d)
    labels = list with measures names

    newLayer = QgsVectorLayer(sourceGeometryType+'?crs='+sourceCRS, name, "memory")
    provider = newLayer.dataProvider()
    attr = sourceLayer.dataProvider().fields().toList()
    attr.extend([QgsField(label, QVariant.Double) for label in labels])
    provider.addAttributes(attr)
    newLayer.updateFields()
    provider.addFeatures(sourceFeats)

After creating the new fields with attr.extend and updating the layer, the dataProvider() method with changeAttributeValues() function is used to populate them with results:
for idxfeat, feat in enumerate(newLayer.getFeatures()):
    featid = int(feat.id())
    for idxlabel, label in enumerate(labels):
        idxfield = int(provider.fieldNameMap()[label])
        val = float(data[idxfeat, idxlabel])
        provider.changeAttributeValues({featid : {idxfield : val}})

QgsMapLayerRegistry.instance().addMapLayer(newLayer)

I'm pretty sure the input passed to changeAttributeValues is correct and the problem is related to the new fields setting. Looks like the dataProvider() can't locate the new fields. I tested on QGIS console fixing the column by one "original" field and it was updated, but trying to use the new created field results in NULL. It doesn't show any error but I could find log messages with what seems to be QGIS DB operations:
Attribute index 17 out of bounds [0;17]
Attribute index 18 out of bounds [0;17]
...

I tested several different methods with layer.startEditing()/commitChanges() but nothing works and people mess up layer edit with dataprovider, this doesn't seems to be a good idea as point here https://gis.stackexchange.com/a/147399/92864. I also checked previous questions like Got Null values when updating attributes with PyQgis and Cannot get changeAttributeValues() to work in a loop but none of them are using memory layer or a copy.
Can someone help me understand what is happening?
Update:
I hardcoded the labels = list with measures names with names copied from the source layer instead of using the new fields and the function updated the columns successfully, confirming that the problem is with the new fields setting. Do I need to set a geometry to the new fields? Adding a field is not just extending a current feature?


texts: Merging rasters which overlap using ArcGIS Desktop?

texts:
I need to merge two rasters which overlap, both indicate distance values. I need the output of the merge to result in the shorter/lower distance value.

I am using ArcGIS 10.1.

I work for a security company. What I am trying to achieve is a map of a building with cameras' resolution. Therefore this map must show the cameras and their field of view so that I can know the resolution in pixel/meter of the camera image at every point of the building. My questions are the following:

I need the distance from every camera to be linear as shown in the pic, not circular as the Euclidean Distance command does as default
I need to keep these output distance rasters limited to the camera field of view, not all around the points
I need to cut in the right way every field of view which overlaps each other (see image)

I'm starting to figure out many things.
I cannot make the ED on all points at once since. Some part of the map can be covered by only camera A but being closer to camera B which is looking at somewhere else; in this case, the ED would give me back the distance from camera B which I am not interested in. I tried to do Mosaic to New Raster with minimum as the mosaic_metod but it doesn't seems to work very well.


texts: Copying annotation feature and pasting within same layer?

texts:
I am having trouble executing a very basic command which is to copy and paste an annotation within the same layer in ArcGIS 10.1.

First, I created an annotation class, then created one label (a balloon call-out with a number) whose attributes I want to copy for all of my locations.

When I went to Editor > Edit Annotation Tool > selected my annotation label > then under Attributes (in my Editor toolbar) > right-click my label > Copy Attributes > then right-click on my Layer > Paste Attributes, nothing happens.

Can someone explain why I cannot paste attributes in this way?