texts: Spatial join many polygon layers and save ones with condition

texts:
I've got 8 polygon layers made from buffering point layers. All have index and geometry columns:

I would like to select areas where at least 3 consecutive buffers intersect (like: buffer_03 & buffer_04 & buffer_05 & buffer_06).

I've got overall over 1 million points so trying to figure an effective way to process this (maybe Python code) as QGIS crashes while trying to Union the layers.


texts: Estimating a Rotated Pole given a gridded dataset with lat/lons at each point

texts:
I have netcdf dataset that I think is regular gridded dataset using a Rotated pole mapping, but the rotated pole information is wrong.

For those not familiar with Rotated Poles, the cartopy documentation has a description of the projection. Check out the source code for the image for a clearer understanding - the Rotated North pole is at 37.5N, 175E (in standard lat/lon coorinates).

The dataset has lat/lon value per grid point, but rotating those into the Rotated pole using the pole data provided results in a nearly regular grid, but it is slightly curved, and rotated. This makes me think that the pole coordinates are wrong.

Is there any simple way to estimate a rotated north pole, given a grid of lat/lons? I guess it should be possible using a linear regression in spherical space over each column of the grid, and then finding the intersection point of all the grid lines, but that seems complex and I'm wondering if there are any pre-existing algorithms for this.


texts: Using Mapbox (LeafletJS), how do I edit imported shapes from GeoJSON?

texts:
Overview: I'm currently working on a project where I'm importing GeoJSON into a map using Mapbox. However, I'm wanting to be able to edit those shapes once imported into the map. I've gotten the GeoJSON imported successfully. I can see my shapes as layers in my FeatureGroup. I can see my shapes on the map as they should be rendered.

The Issue: The problem that I'm having is editing the layers after I've imported them. When I click on the "edit" button using Leaflet Draw, I'm getting an error saying "Uncaught TypeError: Cannot read property 'enable' of undefined" on line 10 of leaflet_draw.js. Also, my shapes are now empty shapes with dashed border and no edit controls. Therefore, I can't even interact with the shapes at all.

My Code: Here's the code that I've used to pull my GeoJSON into the map:

var featureGroup = L.featureGroup().addTo(map);
if (oldData) {
    var featureLayer = L.mapbox.featureLayer(oldData, {
        pointToLayer: function (feature, latlon) {
            if (feature.properties.type === 'circle') {
                return L.circle(latlon, feature.properties.radius);
            } else if (feature.properties.type === 'marker') {
                return L.marker(latlon);
            }
        },
        setStyle: function (feature) {
            return {color: feature.properties.color};
        },
        onEachFeature: function (feature, layer) {
            layer.bindPopup(feature.properties.description);
            layer.openPopup();
        }
    }).addTo(featureGroup);
}

var featureLayer = L.mapbox.featureLayer().addTo(map);
var drawControl = new L.Control.Draw({
    edit: {
        featureGroup: featureGroup
    },
    draw: {
        polygon: true,
        polyline: true,
        rectangle: true,
        circle: true,
        marker: true
    }
}).addTo(map);

Also, it may be helpful to know that I'm using Leaflet.draw version 0.2.3 and Leaflet 0.7.2 that come with Mapbox 2.2.1.

Any tips will help.


texts: A* pgRouting issue with x1,y1 and x2,y2 values

texts:
I am trying to use the A* algorithm in a topology I have made in PostgreSQL. I can fill the x1,y1 and x2,y2 the latitude/longitude, but I have read that the same version has problems with ST_StartPoint & ST_EndPoint and use another functions.

I try but the x1 x2 y1 y2 is stand empty after the execution of my query and give message of success

here are the queries (first I use these four)

UPDATE roads SET x1 = ST_x(ST_startpoint(geom));
UPDATE roads SET y1 = ST_y(ST_startpoint(geom));
UPDATE roads SET x2 = ST_x(ST_endpoint(geom));
UPDATE roads SET y2 = ST_y(ST_endpoint(geom));

after I try these four and the same, nothing happens!

UPDATE roads SET x1 = ST_x(ST_PointN(geom, 1));
UPDATE roads SET y1 = ST_y(ST_PointN(geom, 1));
UPDATE roads SET x2 = ST_x(ST_PointN(geom, ST_NumPoints(geom)));
UPDATE roads SET y2 = ST_y(ST_PointN(geom, ST_NumPoints(geom)));

PostgreSQL does not give any error, but puts an empty column.

As an example after the queries, I have this result

> gid|dromoi_id|type|leng|geom(MultiLineString,4326)|source|target|reverse_cost|x1|y1|x2|y2
> 1 | 10      |as  |300 |     01051232             | 15   | 60   | 352.2       |  |  |  |


texts: QGIS, NAS and "Xerces XML parsing library"

texts:
Since the beginning of this year 2017 the standard geo information and data of Northrhine Westphalia are available for free: https://www.opengeodata.nrw.de/produkte/geobasis/

Including also data from the land registry (land parcels, etc):
https://www.opengeodata.nrw.de/produkte/geobasis/lika/alkis_sek/bda_oe/

I have downloaded and unzipped the data for Solingen and Aachen:
bda_oe_05122000_Solingen_EPSG25832_NAS.zip
bda_oe_05334000_Städteregion_Aachen_EPSG25832_NAS.zip

I use qgis 2.14 - I can select and see the XML-files when I select NAS as data type, which I want to open as vector layer..

But I do not see any content. The attribute table seems to be empty.

I have found out that it needs the "Xerces XML parsing library". But I do not know if it is part of QGIS 2.14?

The is a similar discussion:

Import of NAS-data, (norm-based data exchange interface)

How can I find out if the "Xerces XML parsing library" is installed and if not - how can I install it?


texts: Generating map sheets and index page in QGIS?

texts:
I have a project in QGIS and would like to print it out in multiple pages. I know how to set the pages one by one (pages 1,2,3), but I would like to know if there is an easy way to do the index page with the overview of the whole project?

Here is a picture to clarify my problem:

This overview could also follow in each pages legend and for example highlight the current page.

I'm new to QGIS.


texts:
I found out that it is possible to get the elevation values for particular location, by using Google maps API:

  https://developers.google.com/maps/pricing-and-plans/#details

What I couldn't find is the resolution of these elevation grid values?

Does Google maps use 3 arc-second (90 meters) for any other location outside of USA?


texts: Calculating destination point using sf package

texts:
I have the following data:

lat <- c(-1.404361, -1.403221,-1.387369)
long <- c(35.35373, 35.35047,35.35171 )
distance <- c(220,150,300)
bearing <- c(67,27,159)

df <- data.frame(lat,long,distance,bearing)

I transform and project into an sf object as follows:

df %>% st_as_sf(coords=c("long", "lat")) %>% st_set_crs(4326) -> df

My question is, how do I calculate the destination point for each long/lat origin coordinate, given the distance and bearing? I have used the destPoint in the Geosphere package, but I cannot find the equivalent in sf.


texts: Attaching point grid data to non-grid points

texts:
I have a dataset with loads of fields of environmental (explanatory) data, as a regular grid of points (black dots in image), and another dataset of sampled (response) values as non-gridded points (red dots).

I'd like to append all of the environmental data fields to the sample dataset, but am not sure how. Most of the environmental fields are simple numeric data e.g. depth, temperature, salinity, etc., and subsequently I'm hoping it's possible to have QGis calculate the values at each red dot based on an interpolation from the nearest 'n' black dot sites.

Firstly, if that is possible, then how I can I do it?

Secondly, some of the fields aren't so easy: one is substrate, so the values are character strings e.g. "sand", "gravel" etc, and thus I'd like the red dot field value to take the value of the nearest black dot site.

Four of the fields are compass directions, i.e. numerical from 0 to 360, which won't pose a problem in most cases, but will do for bearings near 360/0, as a two point average of 359 & 0 will equal 179.5 rather than 359.5. Again, using nearest black dot site will very probably suffice.


texts: How to add multiple shapefiles with C# ArcObjects?

texts:
I have a folder contain 5 shapefiles.

I want to develop C# code for adding the 5 shapefiles with single button click.

I have created code for adding a single shapefile but how can I add multiple shapefiles?

ESRI.ArcGIS.ArcMapUI.IMxDocument mxd = ArcMap.Application.Document as ESRI.ArcGIS.ArcMapUI.IMxDocument;
ESRI.ArcGIS.Carto.IMap map = mxd.FocusMap;
ESRI.ArcGIS.Geodatabase.IWorkspaceFactory wsf = new ESRI.ArcGIS.DataSourcesFile.ShapefileWorkspaceFactory();
ESRI.ArcGIS.Geodatabase.IWorkspace ws = wsf.OpenFromFile(@"D:\ARC OBJECTS\Kothagudem Attribute Data", ArcMap.Application.hWnd);
ESRI.ArcGIS.Geodatabase.IFeatureWorkspace fws = ws as ESRI.ArcGIS.Geodatabase.IFeatureWorkspace;
ESRI.ArcGIS.Geodatabase.IFeatureClass fcls = fws.OpenFeatureClass("Agriculture.shp");
ESRI.ArcGIS.Carto.IFeatureLayer lay = new ESRI.ArcGIS.Carto.FeatureLayerClass();
lay.FeatureClass = fcls;
lay.Name = fcls.AliasName;
map.AddLayer(lay);


texts: Dropped geometry_columns and geography_columns tables

texts:
(cross-posted from SO https://stackoverflow.com/q/55869582/2802810)

I've accidentally dropped geometry_columns and geography_columns tables from an existing PostgreSQL PostGIS database schema - let's call it mydb.myschema.

This post outlines how to restore using CREATE TABLE ... with .sql from PostGIS install (I am on a shared server and am having to get my devops to look for them because of permissions).

I wanted to check if this will work. Also, do these two tables (geometry_columns and geography_columns) contain data relating to the geometries held in existing tables in mydb.myschema. Or do geometry_columns and geography_columns hold reference information for spatial operations?

Any advice from anyone who has gone through this before?

Running PostGIS 2.5.2 on PG 11.x.

[NOTE: No backup available as snapshots weren't set up beforehand...]


texts: How to buffer a point in a Python Add-In by ArcMap's selection tolerance in meters?

texts:
As a follow-on to the answer provided by blah238 to this question.

So I'm working in maps with units of meters, and I convert 96 DPI as follows, but I don't seem to be able to click on anything when zoomed beyond a certain scale?

def get_srch_distance_metres(scale, selection_tolerance=3):
    """ returns the map distance in metres that corresponds to the input
        selection tolerance in pixels (default 3) at the specified map scale"""
    # DPI used by ArcMap when reporting scale (96.0 DPI = 3779. Dots Per Metre)
    return scale * selection_tolerance / 3779.527559055


texts:
The contour map that I generate from a Maps Made Easy (MME) DEM produces a distorted result and I’d like to learn how to correct the problem.

Shown is a pool that has been physically measured to be rectangular. To visualize that the edges are orthogonal I’ve rotated the clip and place a rectangle over it (lower left).

Next, I imported the MME DEM into QGIS and generate a 5 ft interval contour map (upper right). Notice that the rectangular pool contour is a rhombus. If the result is exported as a DXF, the DXF file contains a rhombus. Oddly the MME GeoTIFF looks like a rhombus too in QGIS, but not in a TIFF viewer.

To illustrate the distortion, the QGIS contour is rotated and a rectangle has been placed over it (lower right).

I’ve changed the CRS to Google Pseudo Mercator, as recommended by MME, but the distortion does not change.

How can I eliminate the distortion?

Below is the contour line superimposed over the geoTIFF. Even the geoTIFF is distorted.


texts: What are SRIDs like 30008 in an SDE geodatabase?

texts: When I look at the SDE.ST_GEOMETRY_COLUMNS table of my SDE database I'm seeing mostly recognizable four-digit values, but a number of columns have SRIDs like 300006, 300008, and 300047. I can't find descriptions of these anywhere and they aren't recognized when I use an open-source library like Proj4 to transform geometries. Does anyone know what these signify?


texts: How to stop/disable an openlayers 3 animation?

texts:
I have the issue of overlap 2 animations if the previous is really long.

I'd like to stop the current animation for set a new animation with different parameters.

var anima = ol.animation.pan({
        duration: 5000,
        source: map.getView().getCenter()
    });
    map.beforeRender(anima);
}
map.getView().setCenter(ol.proj.transform([lng, lat], 'EPSG:4326', 'EPSG:3857'));

How can I stop the current ol.animation with openlayers 3? I think with duration = 0, but I don't know to access that property :$


texts: QGIS joining two almost identical line features

texts:
I am working on a project in QGIS georeferencing and merging some raster maps that include road data. In the original raster the roads have data assigned to them according to their color in the raster.
I have georeferenced the rasters, used GRASS's r.thin to thin the resulting road lines, and used r.to.vect to convert them to vectors. What I end up with is a set of road lines that are jagged because of the raster pixels.
I would like to take the attributes of this layer (based on the original color) and join them with a more well established and accurate road shapefile that I downloaded from the internet that would add road names, county information, etc.
The picture below shows the two layers. The distance varies across the area of interest but at max is about 40 meters off. The red is the "actual" road shapefile and the black is the raster to vector conversion.
I could use the join attributes by nearest tool, however this joins some raster road attributes to roads that are not shown in the original raster (i.e. neighborhood roads or dead end roads) (see picture). Is there a way to match up these almost identical layers or to join the line features based on an average distance across their length rather than intersections?
On the other hand, would it be easier to go the opposite direction and assign data to the "actual" road vector features by the raster colors they overlap? (see second image) Could I still exclude those dead end roads that are not included in the original raster?


texts:
I need to create a new raster layer with cumulative values for 16 raster layers.  For context, each raster layer represents an 'areas burnt' in a given year - essentially, I want a new layer that tells me how many times each area has been burnt using a colour e.g. white is never burnt, light yellow is burnt once, darker yellow burnt twice… etc.)

I have been told I can use the raster calculator to do this, however when I add the layers together I am (unsurprisingly) left with a raster that only covers where each raster layer overlaps.  Is there a way to get around this?

I actually also have a polygon layer for each raster layer, in case someone has a solution using the polygons.  I converted the polygons to raster because the polygons are quite complicated, and I was told this would make it tricky.  Screenshot of a section of the polygons below for context:


texts: Cut a picture with the grip of another image

texts:
I have two raster images.
Let "A" and "B" my images
with "A" larger than "B"
I would like to cut A for it to have the same grip as B with the python script for use in the automation of a processing chain.
I use the SCRIPT QGIS 2.14 editor. and the modeler QGIS 2.14
I am available to answer all your questions so that you understand better if this is not yet the case.


texts: ArcPad 10.2 Command bar & Quick Access Menu

texts:
I am interested in making the Quick Access menu and Command Bar not available/visible for my users. The Quick Access Menu is prompted by the little icon on the bottom right of the tab bar seen here:

The Command Bar is only visible while making edits but I have the tools it provides already in a custom draw tool bar and I believe it to be unnecessary. Also, it provides a menu icon much like the Quick Access Menu that I do not want end users messing with. For reference the bar is as such:


texts: Opening File GDB in Samba throws exception - but only in VS2010 debug mode

texts:
I'm accessing a file geodatabase through a mapped network drive pointing to a unix machine running Samba.

When I start the app from within VS2010, I get a COM exception when I attempt to open the file Gdb.  No exception is thrown when not running in the debugger.

I'm using the recommended method of opening the gdb (via Activator.CreateInstance):

public static IWorkspace FileGdbWorkspaceFromPropertySet(String database)
{
   IPropertySet propertySet = new PropertySetClass();
   propertySet.SetProperty("DATABASE", database);
  Type factoryType = Type.GetTypeFromProgID(
    "esriDataSourcesGDB.FileGDBWorkspaceFactory");
  IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)
    Activator.CreateInstance(factoryType);
  return workspaceFactory.Open(propertySet, 0);
}

If I copy the same file gdb to the local windows machine, I can access it without any exception.

Has anyone ever run into this issue?

I know very little about Samba - is there something that needs to be configured to make it more file gdb friendly?


texts: Get a column from database and populate the shapefile attributes?

texts:
I would like to know how to get copy a specific columns from the database and add that column in a shapefile attribute table?
I am open to other better ideas.


texts: Create voronoi diagram from line segments

texts:
I'm looking for a way to create a voronoi diagram based on line segments instead of points. See the example below (taken from this presentation).

Ideally I'd like something I can script from Python, but a solution using ArcGIS or similar would also be acceptable.

The only library to do this I've found so far is openvoronoi, which looks promising. Are there any others?


texts: Converting lyrx file or style manager stylx file from ArcGIS Pro to ArcMap

texts:
Is there a way to either convert a lyrx file or the style manager stylx file from ArcGIS Pro 2.0 into anything that can be used in ArcMap 10.5?

I have symbology for about 100 items in ArcGIS Pro that I need to be able to use in ArcMap.


texts:
I have bunch of data points with latitude and longitude. I want to use R to cluster them based on their distance.

I have already taken a look at this page and tried clustTool package. But I am not sure if clust function in clustTool considers data points (lat,lon) as spatial data and uses the appropriate formula to calculate distance between them.

I mean I cannot see how they differentiate between spatial data and ordinal data. I believe the distance calculation between two points on map (spatial) and two normal numbers is different. (Is it not?)

Also what happens if I want to consider a third parameter in my clustering?

Like say if I have (lat,lon) and one other parameter.

How is the distance calculated?

The other problem I have with clustTool is that it is designed with a GUI in mind. I don't know how I can skip the GUI overhead in the library because I don't need it.

What options do I have in R for cluster analysis of spatial data?


texts: Classifying quantiles in raster dataset with ArcPy?

texts:
I am trying to write a Python script which uses the reclassify function in ArcPy to create 9 classes in a raster file. I want to use the quantile method to create the breaks, but so far, I can't seem to find a function that can create quantile breaks in the raster data.

Essentially, I want to perform the above tool with ArcPy. Currently, I have the below in my code, but I do not know how to calculate the remap variable. It cannot be hardcoded as this function should work for all raster inputs. Also, the only modules I have installed are ArcPy and NumPy. I am not able to install any other modules because this is a school project.

tool_input = in_file
reclass_field = "Value"
remap = ???
missing_values = "NODATA"
tool_output = out_files

Reclassify(tool_input, reclass_field, remap, missing_values).save(tool_output)


texts: NWS Grib Files scheduled download script

texts:
I am downloading grib2 files from the NWS NOMADS GFS service details here

The files are provided on a 4x per day cycle and the date format is YYYYMMDDHH. I believe that the processing is actually done 4 times per day.

I am a little confused about the naming convention and what time to actually go and get the files.

An example url using the filter service is:
    http://nomads.ncep.noaa.gov/cgi-bin/filter_gfs.pl?file=gfs.t00z.pgrbf00.grib2&lev_surface=on&var_PRES=on&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%%2Fgfs.YYYYMMDDHH

My understanding is that he HH can be 00, 06, 12, 18.

How does this file naming system work?
Does any one have a script that downloads this data on a scheduled basis and figures out what the file names are?
does this file have forecast weather available in advance?


texts: How to check if a leaflet layer is already added on the map?

texts: I built a simple map app where I have a layer on top of a map. In Leaflet, how can I check if the layer is already present?


texts:
I am supposed to print the area closest to the blue line in the centre of a map that looks like this:

I am not interested in seeing the imagery that is far away from the blue line in the centre. However, it is only supposed to have a scale of 1:10,000, which looks like this in Layout View:

Is there a way to only print the imagery closest to the blue line, but not have it be zoomed in on only one area of the line? I've tried using larger paper but it looks the same on each size.


texts: Save shapefile as grouped/layered vector file

texts:
I want to get the 10m natural earth countries file into a Mercator projection and then into a vector file (SVG, ai, etc) where for each country, its main land and any of its islands are in a single group/layer.

Using QGIS I can change the projection and export a dfx, then open the dfx in Illustrator. But the layers are not grouped by country (e.g. Greece is like 20 different layers - which is fine, but they are not together in a group or layer so there is no way to select "Greece" as a whole).

In QGIS I found the "Select Features by area or single click" button. If I click on Greece using that pointer, it correctly selects ALL of the Greek islands as a group. I see I can export that as it's own file. But I don't want to have to do that 247 times and then resize and puzzle back every country into the world again.

I found Kartograph.py and spent hours trying to install it using the terminal, which I'm not good at, so failed there.

Need some help/suggestions!


texts: Why am I getting two different outputs from Python and C++ using GDAL?

texts:
I am trying to gather information from a *.tif file using the GDAL libraries. When I open the *.tif file using both the Python API, and C++ API, I encouter two completely different answers. Here are the two blocks of code, and their terminal outputs:

Python

import gdal
from gdalconst import GA_ReadOnly

def main():
    data_set = gdal.Open("file_name.tif", GA_ReadOnly)
    print data_set.GetProjectionRef()
    print 'Driver: ', data_set.GetDriver().LongName

if __name__ == '__main__':
    main()

Output

PROJCS["unnamed",PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

Driver:  GeoTIFF

C++

#include "gdal_priv.h"

int main() {
  GDALDataset  *poDataset;

  GDALAllRegister();

  poDataset = (GDALDataset *) GDALOpen("file_name.tif", GA_ReadOnly);

  if(poDataset == NULL)
    return 0;

  if(poDataset->GetProjectionRef() != NULL)
    printf("Projection is `%s'\n", poDataset->GetProjectionRef());

  printf(
    "Driver: %s/%s\n", poDataset->GetDriver()->GetDescription(),
    poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));

  return 0;
}

Output

Projection is

`PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32613"]]'
Driver: GTiff/GeoTIFF

If it helps in the answer. I used the terminal to install the C libraries, and I am using Anaconda to manage my Python libraries.


texts: Problems with gdal1.8 Java bindings gdal.ReprojectImage, produces no data

texts:
I have a problem with gdal.ReprojectImage from gdal1.8 Java
bindings. When I perform the reprojection, the projected image is
totally black. And the origin of the projected image is wrong too.
There were no errors popping up at all. The code is the following:

           Dataset in_ds = gdal.Open(inPath, gdalconst.GA_ReadOnly);
           Dataset out_ds = in_ds.GetDriver().Create(outPath,in_ds.getRasterXSize(), in_ds.getRasterYSize(), in_ds.getRasterCount());

           try {
                   SpatialReference dstRef = new SpatialReference("");
                   dstRef.ImportFromEPSG(26919);

                   out_ds.SetProjection(dstRef.ExportToWkt());
                   out_ds.SetGeoTransform(in_ds.GetGeoTransform());

                   if (gdal.ReprojectImage(in_ds, out_ds) == gdalconst.CE_Failure)
                           System.out.println("something is wrong");
           } finally {
                   in_ds.delete();
                   out_ds.delete();
           }

Does anybody know what's going on here? Thanks a lot.


texts: Installing QGIS with GRASS plugin on Ubuntu?

texts:
I already tried to install QGIS 2.8.3 using the procedure  described on:
http://www.qgis.org/en/site/forusers/alldownloads.html#debian-ubunt,  but grass-plugin do not compare on plugin list!

I tried both repositories: ubuntugis or debian, but nothing!
I also cleaned the previous QGIS 2.6 installation (remove old ppa) on my UBUNTU 14.04.

Does anyone tell me a complete and verified fresh procedure (step by step) to do a good install that allows grass-plugin (grass 6 or 7 it's the same)?

19 ago update :
I've done exactly this steps (before I cleaned previos install and repositories):
    sudo apt-get install python-software-properties
    sudo add-apt-repository 'deb http://qgis.org/debian trusty main'
    gpg --keyserver keyserver.ubuntu.com --recv DD45F6C3
    gpg --export --armor DD45F6C3 | sudo apt-key add -
    sudo apt-get update
    sudo apt-get install qgis python-qgis qgis-plugin-grass grass-gui grass-doc

So, I correctly installed : QGIS 2.10.1, GRASS 6.4.2 and gras-plugin working !!
But ... mysteriously in grass-toolbar:  "add raster layer"  and "add vector layer" do not compare !! ....I d'nt understand.
What's the matter ?

Using repository 'http://qgis.org/debian trusty main', it is installed the Qgis 2.10 release ... if I am not mistaken !

If it's possible instead, I want to install Qgis 2.8.3 LTR with grass-plugin working  (GRASS 6 or 7 it'snt important).


texts: Trying to interpolate depth measurements of a pond

texts: We have taken depth measurements across a pond, and have a CSV file with each location's GPS coordinates and depth in meters. We are trying to create a map of this pond showing where the channel is and what areas have been filled in with silt (shallower). When using QGIS interpolation, the TIN method doesn't seem to work at all, and with the IDW method it is close to what we want, but is very patchy and incomplete. We want a smoother, more complete depth map of this pond. Can anyone explain what is happening in the TIN method picture, or explain how to obtain a better result?


texts: Why do I get different results from PyQGIS feature.geometry().type()?

texts:
On the python console in QGIS what is the difference between the following two code snippets (activeLayer is a polygon layer):

# 1)
iface.activeLayer().selectedFeatures()[0].geometry().type()
>>> 3

# 2)
f = iface.activeLayer().selectedFeatures()[0]
f.geometry().type()
>>> 2

In case 1) the result is '3', in case 2) the result is '2',
and moreover in QGIS 2.12 using 'asPolygon()' instead of 'type()' case 2) gives as list of qPoints as expected while in case 1) QGIS says goodbye with a 'Runtime Error! ... R6025 - pure virtual function call' (in QGIS 2.8 results in an empty list).

I expected the same results.


texts: How to clip ITN Network Dataset in ArcGIS

texts:
I have an Integrated Transport Network dataset that I downloaded from Edina Digimap that I was planning on using for Service Area and Closest Facility analysis in ArcGIS 10.2.2.

I have downloaded, converted and prepared the data for use in ArcGIS following these instructions and the data loads up and can be used for routing etc.

However, I need to clip the ITN Network Dataset so it matches a county boundary that I have as a shapefile so as to restrict routing to within the county only.

The Clip tool won't let me Clip the Network Dataset that was built from the ITN files. I can Clip the individual files themselves (...) but what do I do with the newly clippped shapefiles?

How do I combine the clipped shapefiles to build a new clipped Network Dataset? Or is there a better way to create a clipped network dataset?

I'm lost!


texts: How to do intersection/union across all tables in schema?

texts:
I have a schema with multiple tables, all of which have a particular field, let's call it foo. They all have geometry wkb_geometry and unique key ogc_fid.

I need to create two new tables that would have intersections and unions of all features 'foo' equal to 'bar'.

I can do it for two tables specified explicitly:

-- intersection
create table intersection_results as
  select T1.foo, st_intersection(T1.wkb_geometry, T2.wkb_geometry) as geom
    from T1, T2
     where st_intersects(T1.wkb_geometry, T2.wkb_geometry) and T1.foo = 'Bar' ;

-- union
create table union_results as
  SELECT st_union(wkb_geometry)
    FROM (
      SELECT wkb_geometry FROM T1 where foo = 'bar'
        UNION SELECT wkb_geometry FROM T2 where foo = 'bar'
          ) as moo ;

But what I am struggling with is to create a joined results table across all the tables  in my schema.


texts: Ratio policy for Intersect in ArcGIS Online?

texts:
I have two polygon layers in ArcGIS Online (for organizations) stored as Feature Layers; one is population block, another is the AOI polygon. They are overlapping and I want to find out what portion of the population is "located" within the AOI polygon.

I used to do that in ArcGIS Desktop by creating a Feature Layer with Use Ratio Policy option with Make Feature Layer GP tool and then running Intersect GP tool (which will respect the field ratio settings). An easy thing.

I cannot find a way to do this in ArcGIS Online though. Overlay with Intersect method doesn't split the attribute value; Summarize Within doesn't clip the area and of course no value splitting either. Which tool/workflow should I use?


texts: How to create model maker NDBI indice using Landsat 8?

texts:
I'm having trouble making the model for NDBI (custom made in ERDAS 9.2) using Landsat 8 image?

As a result of running the model I'm getting a white image.


texts: Wrong option for GRASS "v.surf.idw" using PyQGIS

texts:
Using Qgis 2.14.12 and PyQGIS, I get an error message when executing the following:

processing.runalg('grass7:v.surf.idw', layer, 18, 2.0, 'p', True, "%f, %f, %f ,%f"%(1.635, 1.813, 50.896,  51.0205), 0.005, -1, 0.001, "e://out.tif")

The output from

layer.isValid()

is True, but I get this error (warning) in the Python Console :

C:\PROGRA~1\QGIS2~1.14\bin>g.region n= 51.020500 s= 50.896000 e= 1.813000 w=1.635000 res=0.005
   Sorry, <51.020500> isn't a valid option
   Sorry, <50.896000> isn't a valid option
   Sorry, <1.813000> isn't a valid option

This results in a 1 km raster at lat/lon 0/0.

However, when I launch the function from the QGIS toolbox, the output is fine.

I've tried both with GRASS 6 & 7.

Any ideas?


texts: Where is the File Geodatabase Raster Dataset file stored?

texts:
I have a esri file geodatabase with a raster dataset. How can I find out which physical file the dataset is referring to?

(I am retrieving the file geodatabase through the internet. It all seems to go well, but the raster dataset is missing in the file database after retrieval. Zipping the file geodatabase before "transport" solves this problem. To debug the issue I want to see which physical files are missing/different between the zipped and not zipped databases.)


texts: Getting different values for area field calculations using ArcGIS Desktop?

texts:
We have self developed application based on ArcGIS 9.3.1 SP2. Our agents use that application to edit parcels. The parcel layer has a field SHAPE_AREA that is calculated by ArcGIS. Our application gets the value of the field property to show the area of the parcel. However we see now that the area is not always correct. The most flagrant difference we encountered was an area of 17ca (m²) instead of 29ca. I try to find out why the values are not correct and how many of those values are incorrect. To do that I open the shapefile in ArcGIS 10.2.1. and I add a field named "area2" with the same properties as the field SHAPE_AREA (type: double, number format: numeric, precision: 11, scale:4), I use Calculate Geometry to get the area and i use the same projection. Then I add a field named "result" with the same properties. I calculate the field: SHAPE_AREA - area2 to get the differences.

Can anybody tell me why I get different results when both calculations are done by ArcGIS on fields with the same properties? And how we can resolve this problem?


texts: Identify ID of original feature and its surrounding buffer? ArcGIS 10.1

texts:
I have a set of polygon features, stored in geodatabase (ArcGIS 10.1 Desktop).  I want to create a surrounding buffers. I want further to be able to acknowledge the ID of original feature to ID of surrounding buffer. How can I make it?

I though the OBJECTID of original polygon equals OBJECTID of surrounding buffer. That's however not true, as my original data ID does not follow the sequence 0-N, but some numbers are excluded.

How can I make sure to identify the couple: original polygon + its buffer? It there any workaround, ideally in ModelBuilder?


texts:
I am trying to create a contour map from river sonar data. I have converted the .gpx file to a .shp file. I then used the field calculator to change all the soundings, which showed as positive numbers, to negative numbers as it is depth. When I carry out the TIN interpolation. It results with both max and min values showing the same as the mean depth. The
I am using QGIS 3.16. I am pretty new to QGIS but I cant understand what I have done wrong.


texts: Export map as image with selected feature color change QGIS 3.2.0

texts: I changed the color of selected features from the default yellow, to red. However, when I export my map as a png, the color is changed back to yellow. How can I export to an image while keeping the color change of the selected features?


texts: Trouble calculating attribute values in pyqgis

texts:
I recently posted another question about calculating unique IDs for features using python and qgis here Calculating unique ID value for a field in QGIS using python. I am trying to use the same technique of iterating over features from two different layers, check if they intersect, and if they do, calculate the values for a field in one layer using the values from a field in the intersecting layer. The code runs without error, but no values are calculated. It is a string field that I am trying to calculate values for. The b features are polygons and the s features are points. When I step through I find that they are intersecting (as expected) and the NumVal prints as expected. Can anyone help? Here is a snippet of the code:

UPDATE: Based on the answer I have edited the code below, which now functions as expected. I have highlighted the line that missing previously.

                        idxBuffNm = dsslvExplodeLyrWGS84.fieldNameIndex('SiteName')
                        idxCsvNm = shpWnames.fieldNameIndex('SITE_NAME')

                        #Start editing Buffer layer -

this was the missing line:

                        dsslvExplodeLyrWGS84.startEditing()

now follows the rest of the code

                        for b in dsslvExplodeLyrWGS84.getFeatures():
                            for s in shpWnames.getFeatures():
                                if s.geometry().intersects(b.geometry()) == True:
                                    NmVal = s.attributes()[idxCsvNm]
                                    b.setAttribute(idxBuffNm,NmVal)
                                    dsslvExplodeLyrWGS84.updateFeature(b)

                        dsslvExplodeLyrWGS84.commitChanges()


texts: Looking for Floodplain Modelling instructions

texts:
I'm really new to GIS, yet I have a project that I am doing that is GIS entensive.

I want to do a flood plain (for a river whose shapefile i have) and publish whatever output (its a map i presume) on a website for guys to access and consume. Can someone help me with an easy to follow floodplain-modelling instructions?
I would love to follow the tips and learn more.


texts: Why does ManageR and rpy2 not work in QGIS (1.8.0-Lisboa)?

texts:
I have been trying to get R scripts to function and ManageR to work within QGIS (1.8.0-Lisboa). (Windows 7). Python2.7.5. Two problems indicated that R was not working in QGis.

First - when I run an R script there is not output in the pop-up window.
Second - when I click on ManageR (bottom right corner of Qgis) it tells me I should either install rpy2 or check my installation of rpy2.

So, installing (checking rpy2). This in itself has proved very problematic. I found the following posts and 'quick-fix' software:

https://stackoverflow.com/questions/11928277/rpy2-and-python-installation-road-block
https://bitbucket.org/breisfeld/rpy2_w32_fix/issue/1/binary-installer-for-win32

followed the instructions but still to no avail. Sadly I am no computer genius. I can follow simple instructions but am not up to fixing anything myself. So, guys. This is as much a request as a question. (a) how to fix, (b) is anybody working to fix this?

many thanks
BW


texts: Identifying highs on DEM using ArcGIS Spatial Analyst?

texts:
I have a DEM and I want to find each hill and create a polygon for its area, effectively the out closed contour that encloses the hill. I can do this visually but I'm about to work on a larger area so need to find a way of scanning the DEM for structural highs.

Is there a tool that does this, or a workflow that might get me close to my results. I've been looking through the hydrology tools in ArcGIS and they look promising but I'm not sure they can quite produce what I'm looking for - If there was one that would identify areas where lakes would formn that could work if I invert my DEM.

Has anyone done this before of do they know how to do it? I'm using ArcGIS desktop 10.2 at the moment but if you have a solution in other software I'm interested.


texts: Creating point at defined distance along line in PostGIS?

texts:
How do I create a point at a defined distance along a line?

I know we have ST_Segmentize, but it creates more than 1 points - I only need 1.

What I want is:

For each line: Create ONLY 1 point at the first or last x meter (or in a spatial reference system)


texts:
In C++ API if I am not mistaken there is class QgsDefaultValue, so we can programmatically set default values(expression) to a field.

Is there python alternative to that?

All research leads to defining trigger on create feature, but it is not really what I want.


texts: Cannot find plugin for GDAL raster miscellaneous tools ( merge ) in newly installed QGIS

texts:
I am using QGIS 3.8.3.

I'm looking for the Gdal tools "miscellaneous" under the raster menu. It does not appear, nor can I find where to install it in the plugins.

Specifically I am looking for the merge tool within miscellaneous.

Screenshot below. Any help?


texts: Extracting variables from the WorldClim database using raster in R

texts:
I am extracting temperature information using the getData function from the raster package.

bio1 <- getData('worldclim', var="bio", res = 0.5, lon=1, lat=-77)

bio2 <- getData('worldclim', var="bio", res = 0.5, lon=7, lat=-70)

The variable bio var=bio , means it is extraction 19 layers of different bioclimatic variables at the lon/lat coordinates provided.

The problem is that if I get the mean of these layers for both bio1 and bio2, even though they are from different lat/lon coordinates, they are the same. And I am using the finest resolution res = 0.5 :

> cellStats(bio1,stat = "mean")
   bio1_46    bio2_46    bio3_46    bio4_46    bio5_46    bio6_46    bio7_46    bio8_46    bio9_46   bio10_46   bio11_46
 159.38445  146.23495   54.29901 4262.23057  292.49362   23.63659  268.85703  178.42543  132.43352  211.26331  102.22708
  bio12_46   bio13_46   bio14_46   bio15_46   bio16_46   bio17_46   bio18_46   bio19_46
 416.39054   61.81748   13.82139   45.91476  165.43986   48.51057  129.29831   72.57548
> cellStats(bio2,stat = "mean")
   bio1_46    bio2_46    bio3_46    bio4_46    bio5_46    bio6_46    bio7_46    bio8_46    bio9_46   bio10_46   bio11_46
 159.38445  146.23495   54.29901 4262.23057  292.49362   23.63659  268.85703  178.42543  132.43352  211.26331  102.22708
  bio12_46   bio13_46   bio14_46   bio15_46   bio16_46   bio17_46   bio18_46   bio19_46
 416.39054   61.81748   13.82139   45.91476  165.43986   48.51057  129.29831   72.57548

From the manual I understood that cellStats would get the mean for the cells of each of my layers. But maybe it is doing something different?


texts: Changing parameter values according to other parameter choice in Processing script using QGIS 3?

texts:
I couldn't find anything about this problem.

Suppose to have the first parameters in a Processing script defined as an Enum:

   self.addParameter(
        QgsProcessingParameterEnum(
            self.CHOICE,
            self.tr('Choose the parameter'),
            ['param1', 'param2']
        )
)

While the second one as a number:

self.addParameter(
    QgsProcessingParameterNumber(
        name=self.VALUE,
        description=self.tr('My value'),
        type=QgsProcessingParameterNumber.Double,
        defaultValue=40,
        optional=False,
        minValue=0,
        maxValue=100
    )
)

Is it possible to change the minValue, maxValue, description, ecc.. of the second parameter depending on the first or second choice of the first parameter already in the GUI?

e.g. if the user choose param1 then the second parameter should become:

self.addParameter(
    QgsProcessingParameterNumber(
        name=self.VALUE,
        description=self.tr('My value'),
        type=QgsProcessingParameterNumber.Double,
        defaultValue=40,
        optional=False,
        minValue=0,
        maxValue=100
    )
)

while if he choose param2 it should became:

self.addParameter(
    QgsProcessingParameterNumber(
        name=self.VALUE,
        description=self.tr('My AWESOME value'),
        type=QgsProcessingParameterNumber.Double,
        defaultValue=0,
        optional=False,
        minValue=50,
        maxValue=500
    )
)


texts: unfortunately I have tried to download images bulk MODIS (NDIV and LST), but since the USGS server ((http://e4ftl01.cr.usgs.gov/) exchange ftp http, based codes R and python obsolete because I know are based on fTP, I tried using pymodis but I could not install I followed all the steps of the page pymodis (http://www.pymodis.org/) but still I could not can someone help me install it or a way to download bulk data ??


texts: Creating 4 band image from RGB and grey scale image in QGIS?

texts:
I have an RGB image and a Gray scale image. I would like to merge them in a 4 bands image.

How can I do that?

Is there a simple way to do it from QGIS or I have to modify the image before to load them in QGIS.


texts: Convert shapefiles from WGS84 to Miller Projection

texts: I have some shapefiles of the world from Natural Earth at WGS84.  I would like to convert them to the Miller Projection.  I have ogr2ogr, but I'm not sure what the parameters  for this conversion should be, or if this conversion can be done with it.


texts:
I want to use algorithms from GRASS GIS in R. One essential step is already providing some serious problems. I managed in some way to install GRASS on Ubuntu 18.04, but now I need to provide the path to the initGrass()-function from rgrass7. I tried using the find2GI-package to find my GRASS GIS installation but it gives me this error:

> link2GI::findGRASS()
Error in gregexpr(pattern = "\"", root_dir)[[1]] :
  subscript out of bounds

Is there any easy way to just find the path to the binaries of grass. I tried which grass aswell but it gives me:
/usr/bin/grass, what also doesn't seem to work with initGrass


texts: Running PostGIS ST_ClusterWithin on large dataset?

texts:
I am trying to use the new ST_ClusterWithin in PostGIS 2.2 to help me aggregate building polygons that are within 5 metres. The source data has > 13million polygons.

I ran the following query on a sample and it worked brilliantly.

CREATE TABLE output AS SELECT ST_CollectionExtract(unnest(ST_ClusterWithin(geom, 5)), 3)::geometry(Polygon, 27700) as geom FROM building

however trying to run the same on the full set causes the following error.

ERROR:  array size exceeds the maximum allowed (1073741823)

Hvaing done googling on that error it seems it is causes by trying to keep the results in memory.

Does anyone know how I could improve my query to work on such a big dataset?


texts: Are there any standard, publicly available symbol sets for crime mapping?  What are their intended applications and their strengths and weaknesses?


texts: Marker points don't fit with street (lack of accuracy while zoom)

texts:
I'm using OpenLayers 3.x with a Google Maps layer (cannot give you link, search "ol3 Google Maps layer").

I can't enhance the accuracy how points would be represented on map, see images below:

This is the zoom out for a series of points. Each of them were accurately took from Google map (right click --> "What's here?"). Then I add to a source layer like this:

var marker = new ol.Feature({
        geometry: new ol.geom.Point(ol.proj.transform(myCoord, 'EPSG:4326', 'EPSG:3857')),
        name: _markerInfoEnrichment(featureItem)
    });
source.addFeatures(marker);

If you pay a little attention, you can actually notice that some points (like top-left of map) has not been accurately represented. Indeed they are rerouted away from street.

But in the next image we can see that the more you zoom the more accuracy you gain:

And if I continue to zoom, the point will be well represented on street.

So, there is a way to enhance the points representation accuracy? Please also notice that same issue is still there with OpenStreetMap too.


texts: Determining if feature is editable using ArcObjects?

texts:
I am using the FocusMap.FeatureSelection to return all the currently selected features.

I want to process only the currently editable features and leave the remainder alone.

How can I quickly programmatically tell which features are editable?

I know of the IEditLayers.IsEditable method, but think that might be a little heavy handed to try and map from a feature back to the featurelayer and then test that for all selected features.

Is there another way, or is this all I have available to me?


texts: Google Earth stopped working on Debian testing

texts:
Suddenly GE Pro desktop stopped working on Debian testing laptop. I can launch GE, menus work. I can even search for a place and pop up appears but where a map shall be only black screen is. Also there are no items in Places and Layers.

Here is my GPU info:

inxi -G
Graphics:
  Device-1: Intel Mobile 4 Series Integrated Graphics driver: i915 v: kernel
  Display: x11 server: X.Org 1.20.7 driver: modesetting unloaded: fbdev,vesa
  resolution: 1440x900~60Hz
  OpenGL: renderer: Mesa DRI Mobile Intel GM45 Express v: 2.1 Mesa 19.3.3

I figured out I'm probably using inappropriate GPU driver but I'm not sure how to solve this problem.

So far I tried to reinstall and downgrade GE, but previous versions depend on lsb-core which is not in use in Debian buster. I also used all option in Launch repair tool. No advance.


texts: ArcGIS Feature Service - Return a Single Row / Record

texts:
Is it possible to get just one row/record returned from a feature service query -- similar to the TOP, ROWNUM and LIMIT in database queries? Of course, OBJECTID "1" is not guaranteed. Is there a trick that I've missed?

http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/SanFrancisco/311Incidents/FeatureServer/0/query?where=OBJECTID%3D1

SQL Server / Access: SELECT TOP 1 * FROM CUSTOMERS;
Oracle: SELECT * FROM CUSTOMERS WHERE ROWNUM <= 1;
MySQL: SELECT * FROM CUSTOMERS LIMIT 1;
PostgreSQL: SELECT * FROM CUSTOMERS LIMIT 1;


texts: Spatialite database creation problem (QGIS DB manager)

texts:
I am trying to create a Spatialite database in QGIS 1.8.0. Step one of instructions is:

Right click on the Spatialite entry in the Browser tree and select Create Database.

The Spatialite icon is frozen. It does not respond to right click or anything else.
Anyone who encountered this issue before?


texts:
I've downloaded some 3D data from 3dwarehouse.sketchup.com

Now, I would like to open it in ArcGIS Pro but I cannot find how can to do it.

How can I do this?


texts:
I have a huge data set with geolocation information. I am trying to visualize the data but the tools cannot handle such a big graph. Size of the file is about 5GB and has more than 3 million rows.

What is the best way to visualize the time-series and geolocation based data?

I have a CSV and I did try using online tools like CartoDB, Plotly and did try using ArcGIS but that is commercial tool.

I got stuck as the file is too big and the tools won't allow loading the files. It is a social network data with spatio-temporal information. I used gepi, cytoscape, ScoNet, Tupil, Anaconda GraphViz but nothing would even load the file. The same if I use sample 10 rows it would work good.


texts:
Looking at this map data from naturalearth (part of Oman & UAE to be specific), each country is defined as a multipolygon. In this case UAE (in red) contains two polygons and Oman (beige) three.

Problem is that QGIS claims UAE to be invalid, since the smaller polygon is "inside" the other one. I've tried making the small polygon a ring instead, but that's also not okay, since it does not touch the other ring.

How should this be organized in order to be valid? Must these be separate features? Joined physically?


texts: How to add a Google Maps layer to Openscales?

texts: I'm using Openscales in a web mapping project. I would like to add a Google Maps layer, but I haven't found any snippets or examples. Any solution will be very helpful.


texts:
I have about 1500 .txt files within a folder named 3410001ne => 3809962sw.  I need to subset about 470 of these files to process with a Python script.  Below is the section of code prior to my for loop which lists all of the files in a folder:

# Get the list of the text files to process
txt_list = arcpy.ListFiles("*.txt")

How can I subset the folder contents to include 470 of 1500 files ranging from 3609902sw => 3610032sw?


texts: Soil maps - displaying specific soil types QGIS

texts:
I'm trying to isolate a couple specific soil types from the GIS soils layers for a county in my state in the US.

I downloaded the data:
ftp://ftp.igsb.uiowa.edu/gis_library/Counties/Adair/soil_01/soils_01.zip

And loaded the .shp file as a vector layer. I can see all the shapes of different soil types but can't figure out how to isolate and display just a couple of them.

Can anyone help with this problem?


texts: the layer attribute table is not consistent with the pointvalue

texts: The point value of my grid is right, but the output value ranges from -2147483648 to 2147483647. I haved tried to recalucate statistics, but it didn't work..


texts: Huge differences between results for 7-param datum transform

texts:
I'm trying to transform a WGS84 lat/lon point

11d42'32.10629" E 5d12'56.75371" S

in three different software packages (Proj4, GeoTrans and Leica GeoOffice) but I receive very big differences between the results - around few hundred meters in X and Y! Note that these variations appears only with 7 parameters not with 3.

Proj4

  cs2cs
  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to
  +proj=utm +zone=32 +ellps=clrk80 +towgs84=178.3,316.7,131.5,-5.278,-6.077,-10.9792,-19.166 +south +units=m +no_defs

GeoTrans

Delta X 178.3
Delta Y 316.7
Delta Z 131.5
Rot X -5.278
Rot Y -6.077
Rot Z -10.9792
SF = -19.166 / 0.999980834 (0.000019166)

Leica GeoOffice

The results are respectively:

800392 9422525
800306 9422840
800941 9422891

I believe all the 3 packages use same mathematical methods for 7-param transform (Bursa-Wolf method). What could be the cause of this huge variation?


texts: SEGFAULT when creating own index in postgresql for POINT3DZ

texts:
I'm trying to write my own index for 3d points - POINT3DZ (using pgis and postgres 9.4). I have written a code similar to spgquadtreeproc.c [1]

This is my spg_quad_point3dz_picksplit function:

PG_FUNCTION_INFO_V1(spg_quad_point3dz_picksplit);
/* Decides how to create a new inner tuple over a set of leaf tuples. */
Datum
spg_quad_point3dz_picksplit(PG_FUNCTION_ARGS){
    spgPickSplitIn *in = (spgPickSplitIn *) PG_GETARG_POINTER(0);
    spgPickSplitOut *out = (spgPickSplitOut *) PG_GETARG_POINTER(1);
    int i;
    POINT3DZ centroid, ctemp;

#ifdef USE_MEDIAN
    /* Use the median values of x and y as the centroid point */
    POINT3DZ      **sorted;

    sorted = palloc(sizeof(*sorted) * in->nTuples);
    for (i = 0; i < in->nTuples-1; i++){
        datumGetPoint3dzP(in->datums[i], sorted[i]);
    }

    centroid = palloc(sizeof(*centroid));

    qsort(sorted, in->nTuples, sizeof(*sorted), x_cmp);
    centroid->x = sorted[in->nTuples >> 1]->x;
    qsort(sorted, in->nTuples, sizeof(*sorted), y_cmp);
    centroid->y = sorted[in->nTuples >> 1]->y;
#else
    /* Use the average values of x and y as the centroid point */

    for (i = 0; i < in->nTuples; i++)
    {
        datumGetPoint3dzP(in->datums[i],&ctemp);
        centroid.x += ctemp.x;
        centroid.y += ctemp.y;
    }
    centroid.x /= in->nTuples;
    centroid.y /= in->nTuples;
#endif

    out->hasPrefix = true;
    out->prefixDatum = (Datum)&centroid;
    out->nNodes = 4;
    out->nodeLabels = NULL;     /* we don't need node labels */
    out->mapTuplesToNodes = palloc(sizeof(int) * in->nTuples);
    out->leafTupleDatums = palloc(sizeof(Datum) * in->nTuples);

    for (i = 0; i < in->nTuples; i++)   {
        POINT3DZ p;
        datumGetPoint3dzP(in->datums[i], &p);

        int quadrant = getQuadrant(&centroid,&p);

        out->leafTupleDatums[i] = (Datum)&p;
        out->mapTuplesToNodes[i] = quadrant- 1;
    }

    PG_RETURN_VOID();
}

Constans USE_MEDIAN is not defined so in the function are executed lines from top of function, next - lines from else to endif, and finally rest of code. Unfortunately first for loop work fine, but second take an error on executing datumGetPoint3dzP.

This is how the datumGetPoint3dzP function looks inside:

void datumGetPoint3dzP(Datum dt, POINT3DZ *p){
    GSERIALIZED *geom;
    LWGEOM *lwgeom;
    LWPOINT *point = NULL;
    geom = (GSERIALIZED *)PG_DETOAST_DATUM(dt);
    if ( gserialized_get_type(geom) != POINTTYPE )
        lwerror("Argument to X() must be a point");
    lwgeom = lwgeom_from_gserialized(geom);
    point = lwgeom_as_lwpoint(lwgeom);
    getPoint3dz_p(point->point, 0, p);
}

I found in debugging that an error occurs while calling getPoint3dz_p.

getPoint3dz_p is a function from PostGIS (lwgeom) [2]

ERROR:

CREATE INDEX example_sensors_index1 ON sensors_data USING spgist( geom geometry_point3dz_ops );
DEBUG:  StartTransactionCommand
DEBUG:  StartTransaction
DEBUG:  name: unnamed; blockState:       DEFAULT; state: INPROGR, xid/subid/cid: 0/1/0, nestlvl: 1, children:
DEBUG:  ProcessUtility
DEBUG:  Creating index "example_sensors_index1" on table "sensors_data"
server closed the connection unexpectedly
This probably means the server terminated abnormally
before or while processing the request.

LOGS:
Some of the messages are in Polish. To everyone I posted a translation under the message.
Here is a part of my logs: https://pastebin.com/dYyb0tFb

How do I fix this code?

links:
[1] 37.187.79.115/postgresql-9.4beta3/src/backend/access/spgist/spgquadtreeproc.c
[2] http://postgis.net/docs/doxygen/2.1/da/de7/liblwgeom_8h_acfc410375b25ddd41feec1b7337173c6.html


texts: turf.js point never found inside or within GeoJSON features loop

texts:
I have a GeoJSON http://radikal.com.mx/distrito.json and all I am trying to do is find if a turf point that has a longitude and latitude is inside or within that GeoJSON. No matter what point I use I always get false.

This is based on the work from: Bill Chappell

Here is my Javascript:

// FILE: js/get-data.js
$(function() {
  var cluster = [];
  var features = [];
  var point = turf.point([-98.768159, 20.098287]); // Hgo | This should return one TRUE at least
  // var point = turf.point([-79.601038, 43.6565353]); // Toronto | This should always be FALSE
  var url = 'http://radikal.com.mx/distrito.json';

  $.getJSON(url, function(data)
  {
    features = data.features;

    for (var i = 0, len = features.length; i < len; i++)
    {
      // Get coordinates
      var is_inside = turf.inside(point, features[i]);
      console.log(is_inside);

      if (is_inside) {
        console.log('Found it inside');
      }
    }

    console.log(cluster);
  });
});

Here is my HTML

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width,initial-scale=1">
  <meta http-equiv="X-UA-Compatible" content="IE=edge">
  <title></title>


texts: Getting properties of each feature in feature collection using Google Earth Engine?

texts:
I would simply like to get the properties of each feature in a feature collection. For example, in the screenshot, I would like to return the value of "1" from "site."

Here is my code:

var table = ee.FeatureCollection("users/adrianom/sites_aea");
var SearchArea = table;
var size_fc = SearchArea.size().getInfo();
var fc_list = SearchArea.toList(size_fc);
var first_fc = ee.Feature(fc_list.get(0));
var fc_properties = first_fc.propertyNames();

This will get the property names of a single feature, but I want to loop through every property. Additionally, I would like to know if there is a better way of getting to a single feature's properties in a feature collection.


texts:
I have a world map. I need randomly to generate polygons and render them on map.
I have no problem with rendering, but can anybody suggest how to generate polygons with random geometry? (I need n polygons, n - random number). The main problem is how to generate geometry to them.


texts: Spatial graphing, distribution of polygon area over latitudes

texts: Given a polygon in a PostGIS database, what is the best way to create a graph of its distribution over a range of latitudes? I'm picturing a graph of a curve, where the proportion of the total area under the curve that lies between two points on the x-axis (latitude) represents the proportion of the polygon's area between those latitudes. The tools I have at my disposal are PostGIS and R. It seems like this should be a fairly routine exercise for which there are well established approaches that I simply haven't been able to find.


texts: Multiple Text rows in QGIS divided by lines

texts:
My Problem seems rather simple but I had no sucess for over a day now.
Basically I have multiple (polygon) shape files that need text information in them (or if the shape is to small outside of it with a line pointing to where it belongs) and show in the map. Its only for 10 layers so I wouldnt mind doing it manually - if i knew how. For a similiar task I added invisible point layers and labeled them (i know probably should have used labeling tool) but this problem remains tricky since the text information is rather complex.
The text hast to look something like this.

Sometimes only two rows, sometimes 4, different information in them per layer. I cant find any option to do multiple rows, that are divided by this lines at all.
I dont need a complex solution, If there is a possibilty where I can enter this stuff by hand as text I would be happy. But the format hast to be like this which is my problem.


texts: Changing order layers are displayed in using MMQGIS?

texts: I am dealing with hundreds of different layers, all of which have one of 6 names - not by choice, that's just the way they show up when I import them.  I am trying to merge all those with the same name using MMQGIS>Combine>Merge Layers, but the layers are listed in a random order in the Merge Layers window so I cannot block off the layers I want to select using shift-arrow.  The order the layers are displayed in is different to the order in the layers panel.  I can press control & click to select different layers, but this is very time consuming and if I make a mistake all the previous selections disappear and I'm back to square one.  How can I change the order the layers are displayed in?  Click and drag would take forever.  Alphabetical order would solve my problem.


texts: In ArcGIS Modelbuilder: Mean value of different polygons within a buffered zone, weighted by area percentage

texts:
I am creating a model in ArcGIS model builder (10.2). I have buffered zones (polygons) and within these zones there are different soil types (polygons; Shapefile feature class). I have assigned values (ranging from 0 to 10) to the different soiltypes, depending on their usability for the concrete purpose (usability score).

What I want to do now is to calculate then MEAN usability score within a buffered zone, weighted by the the share of the land.
For example, if a soiltype polygon with usability score 10 covers 100% of the buffererd zone, the MEAN usability score would be 10, but in case it only covers 50% and a soiltype polygon with a usability score of 2 covers the other 50%, the MEAN usability score should be 6.

To sum up, I am looking for something like "Zonal statistics" (which exists for raster layer as input value), but adaptable to polygons with a certain score as value.

I could not find any solutions yet, maybe there is a standard operational procedure I do not know about?


texts: Generating polygons that represent intercardinal quadrants (NE, SE, SW, NW) at Road intersections using ArcMap

texts:
I have a road polyline feature class in an Oracle 18c/10.7.1 geodatabase.

I would like to generate a polygon FC that represents the quadrants of the road intersections.

The polygon FC would have a QUADRANT column.
Each polygon would have a value that represents the approximate intercardinal direction of the polygon (NE, SE, SW, NW).
The road polylines are angled at -19 degrees.

Is there a way to do this using ArcGIS Desktop 10.7.1 (advanced license)?


texts:
I have a triangle and an arbitrary polygon lying on the same plane in 3D space. I need to find the fastest way to determine if these two intersect, one is inside the other, or completely separated. At first stage, it's not important to find the intersection points or tell the difference between the states when one is inside the other and when they intersect.

The important thing is having the minimum number of calculations. Can someone guide me in the right direction?


texts: Using If/Then Statement in ArcGIS Field Calculator?

texts: I am trying to add a column that shows vegetation production (numeric) in relation to soil series. I have each soil series already in the attribute table however I am having a hard time writing the if/then statement to say something along the lines of "If soil series is "vc" Then production should be 2500". Can someone help me?


texts: Catch the VisibilityChanged event for TopologyLayer

texts:
I can set up a handler to catch the VisibilityChanged event of every map layer returned by:

pEnumLayer As IEnumLayer = My.ArcMap.Document.ActiveView.FocusMap.Layers(Nothing, False)

Then with each layer returned this works:

pLayerEvents As ILayerEvents_Event = pMapLayerReturned

Unless the layer is a TopologyLayer. If it is a TopologyLayer it's an invalid cast. It appears that TopologyLayer class doesn't use ILayerEvents. How then can I tell when the visibility of a TopologyLayer changes?


texts: How to load map tiles according to the drag/touch event

texts:
We are building map apis for javascript integrated with our own map data server(tiles). (Though we can use openlayers or something else, but we can not for some thrid-part reasons).

In a word, we are re-writing something like openlayers.

The following is what we have done:

When a map is created, it will have the container(dom element),center(LatLng), and zoom . Then we will orginaize the tiles inside the container according to the center and zoom.
When user drag the container. We will record the increasement of the mouse position. Then caculate the new center for the map, and then re-caculate the required tiles and orginaize them.

Like this:

onMouseMove(dx,dy){
  var dinx=caculateDegreeInCurrentResolution(dx);
  var diny=caculateDegreeInCurrentResolution(dy);
  this.center=new LatLng(this.center.lat+diny,this.center.lnt+dinx);
  this.refresh();
}

It works as expected, but we found that it is hard to drag the map. It is not smooth like openlayers or other map library during the drag operation.

So I wonder if we are going the wrong way?


texts:
I wonder if there is any publically available dataset that has left turn signal information (say, in the U.S.)? For example, in right-hand traffic countries, major intersections may have a left-oriented arrow to protect left-turning cars.

I've looked into the traffic_light tag in OSM, but seems there is no left-turn protection information.


texts: Finding the closest facility in a feature layer from a webmap

texts:
I have created a webmap on which a feature layer is displayed and I'm trying to find the closest facility when the user clicks on the map.

Based on Find closest facilities in the Javascript API reference, I used the data given by map.getLayer().graphics to initialize my facilities set (http://jsfiddle.net/nikoshr/ob5rxsd7/ for a demo, http://bit.ly/14VEwL3 for the map):

arcgisUtils.createMap("9a33727d54904de7b80687fb0810119a", 'map').then(function (response) {
    var map = response.map;

    var lay = map.getLayer(map.graphicsLayerIds[3]);
    var facilities = new FeatureSet();
    facilities.features = lay.graphics;

    var params = new ClosestFacilityParameters();
    params.travelMode = 5;
    params.facilities = facilities;
    params.returnIncidents = true;
    params.returnRoutes = true;
    params.returnDirections = true;
    params.outSpatialReference = map.spatialReference;

    map.on('click', function(evt) {
        var inPoint = new Point(evt.mapPoint.x, evt.mapPoint.y, map.spatialReference);
        var location = new Graphic(inPoint);

        var incidents = new FeatureSet();
        incidents.features = [location];
        params.incidents = incidents;

        closestFacilityTask = new ClosestFacilityTask("https://route.arcgis.com/arcgis/rest/services/World/ClosestFacility/NAServer/ClosestFacility_World");

        closestFacilityTask.solve(params, function(solveResult){
            console.log(solveResult);
        });
    });
});

and I get the a response looking like this

{
    declaredClass: "esri.tasks.ClosestFacilitySolveResult",
    directions: [{routeId:1, routeName:"Location 1 - Point 1", summary:{...}, ...}],
    facilities: null,
    incidents: [{geometry:{...}, attributes:{...}, declaredClass: "esri.Graphic", ...}],
    messages:[
        {type:50, description:"Input field [VISIBLE] was not mapped to a field in the network analysis class Facilities.", declaredClass:"esri.tasks.NAMessage"},
        {type:50, description:"Input field [TYPEID] ...", ...},
        {type:50, description:"Location Point 3 in Facilities is on a soft-restricted network element.", declaredClass:"esri.tasks.NAMessage"}
    ],

    pointBarriers: null,
    polygonBarriers: null,
    polylineBarriers: null,
    routes: [{geometry:{...}, attributes:{...}, declaredClass:"esri.Graphic", ...}]
}

Onto the questions:

The routes and directions attributes seem correct, but how can I get the facility used to build them? I currently use this but it relies on kowing what was sent to the server

var route = solveResult.routes[0],
    fid = route.attributes.FacilityID,
    facility = facilities.features[fid-1];

Is there a way to pass an url representing the feature layer for the facilities  as hinted in the Closest facility service doc, instead of copying them and sending them back as I do now? Something like

var facilities = new FeatureSet();
facilities.url= 'an url describing the features';


texts: Draw GeoJSON polygon over image returned by Mapbox in Python

texts:
I am trying to draw a polygon representing the GeoJSON on top of the image I get from a map returned by Mapbox. If I pass in the GeoJSON in the Mapbox REST API URL, I get the following image.
.
But I want to draw the polygon manually using PIL as there is a size 8kb restriction for Mapbox API (https://docs.mapbox.com/api/maps/static-images/#overlay-options).
If I just get the map and draw the image using the following code, my polygon seems to be in the correct shape but too small. Can someone point me in the right direction on how to get the polygon size correct?
    def draw(self, mapbox_image, zoom, **kwargs) -> None:
      geometries = get_geometries(kwargs['geojson'])
      left_x = kwargs['left_x']
      top_y = kwargs['top_y']
      actual_x_padding = kwargs['actual_x_padding']
      actual_y_padding = kwargs['actual_y_padding']

      draw = ImageDraw.Draw(mapbox_image)
      for geometry in geometries:
         for coordinates in geometry['coordinates']:
            polygon = [( actual_x_padding+ (longitude_to_x(coordinate[0], zoom ) - left_x),
                       (actual_y_padding +latitude_to_y(coordinate[1], zoom ) - top_y)) for coordinate in
            coordinates]
            draw.polygon(polygon, fill=None, outline="blue")

    ##This is how the parameters are calculated before passing in to the function

    zoom = get_required_zoom(extent_rectangle_lat_long, IMAGE_WIDTH, IMAGE_HEIGHT, IMAGE_PADDING)

    left_x = longitude_to_x(extent_rectangle_lat_long[0], zoom)
    right_x = longitude_to_x(extent_rectangle_lat_long[2], zoom)
    top_y = latitude_to_y(extent_rectangle_lat_long[3], zoom)
    bottom_y = latitude_to_y(extent_rectangle_lat_long[1], zoom)

    centre_x = left_x + ((right_x - left_x) / 2)
    centre_y = top_y + ((bottom_y - top_y) / 2)
    centre_longitude = x_to_longitude(centre_x, zoom)
    centre_latitude = y_to_latitude(centre_y, zoom)
    actual_x_padding = 2 * (longitude_to_x(extent_rectangle_lat_long[0], zoom) - image_left_x)
    actual_y_padding = 2 * (latitude_to_y(extent_rectangle_lat_long[3], zoom) - image_top_y)

I use the calculations given in MapBox or Google Maps Static API: Extent or bounding box of image? to for the conversion.
def latitude_to_y(latitude_degrees, zoom):
  pi_by_4 = math.pi / 4
  latitude_radians = math.radians(latitude_degrees)
  return (128 / math.pi) * math.pow(2, zoom) * (
        math.pi - math.log(math.tan(pi_by_4 + (latitude_radians / 2)), math.e))

def longitude_to_x(longitude_degrees, zoom):
    longitude_radians = math.radians(longitude_degrees)
    return (128 / math.pi) * math.pow(2, zoom) * (longitude_radians + math.pi)


texts: How to change tileSize in XYZ from default([256,256]) to something else for openlayers?

texts:
I am trying to change the tile size and I have not been able. Below is the code.

var projection = get('EPSG:4326');
local_url = 'assets/tile/{z}/{x}/{y}.png';

map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new XYZ({
        url: this.local_url,
        projection: projection,
        wrapX: false,
        tileSize: [1024,1024]
      })
  ],
  view: new View({
    center: [0,0],
    zoom:0,
    minZoom: 0,
    maxZoom: 5,
    projection: projection
  })
});

No matter what I set the tileSize to be it is always [256,256].


texts:
I'm getting this output as a result of discharge values while doing flood inundation mapping.

In some areas there are breaks as indicated by red arrows.

If this result is correct what would be the interpretation of it?


texts: Failing to Delete Vector File from Personal Geodatabase?

texts:
I have a personal geodatabase - based on an MS Access .mdb file - from which I cannot delete an old feature class.

The tables relating to the feature class have been deleted from the .mdb, but the feature class remains visible as part of the geodatabase in ArcCatalog 9.3.

I cannot delete it through ArcToolbox, nor through ArcCatalog, either by right-clicking, or selecting the feature class and choosing the 'delete' button.

I receive this error message:

  'Failed to delete selected object(s) The application is not licensed
  to perform this operation The table was not found. The Microsoft Jet
  database engine cannot find the input table or query
  'postcodes_Project'. Make sure it exists and that its name is spelled
  correctly.'

I understand that the program is struggling as a result of the tables being missing from the database, but is there anything I can do to delete it permanently?


texts: Checking for the presence of a particular field in ModelBuilder

texts: The model that I'm now building in ModelBuilder will take as an input a certain feature class that may, or may not, have a certain field. I would like the model to check for the presence of that field and select only certain records by it if it exists. If the input doesn't have that field, the model is to skip this step and proceed straight to what it what do afterwards. At If statement in ModelBuilder to check if field exists, the Python script for something similar is discussed, but I wonder if it is possible to do this in ModelBuilder somewhat differently--i.e., export the selected features to a new feature class if the field exists. What tools should one use for this if this is possible?


texts: Reprojected Shapefile does not align with OpenLayers Plugin

texts:
I am based in NSW, Australia.  I have a shapefile that when I underlay with the OpenLayers plugin (Google satellite), the shapefile is skewed despite all my efforts to correct it.  I changed the project properties to GDA 94 and GDA 94 MGA 56, as well as re-saving the shapefile as WGS 84 EPSG:3857 to match the OpenLayers plugin.  Nothing seems to work.

Can anyone help?  Is it the original shapefile?


texts: Measuring straightness of a curve segment (represented as a polyline)

texts:
I'm working on an automatic elevation contour labeling algorithm and one of the factors that I want to take into account when deciding the positions of labels is how "straight" a particular segment of a contour is. The more straight it is, the more likely it is to be used to place the label on that segment.

Each contour is represented by a polyline (but with points close together as to look like a curve to a naked eye). I then have a fixed length (width of a label), say, 100 pixels. If I randomly (or otherwise) choose a contour segment with the width of 100 pixels, I want to be able to get a numerical quantitative value of its straightness (say zero for a totally straight contour segment, some value larger than zero for a not so straight segment, and this value increasing as the crookedness increases).

I've searched around for answers but I couldn't find anything really useful. I'd be grateful for any pointers.


texts: Determining Upslope and Downslope from Line Features to Map Hydrologic Capture Zones

texts:
I am trying to map hydrologic capture zones (HCZ).  HCZ's are contributing areas upslope of any potential capture feature (i.e. road, pipeline path, railroad line), that would alter/impact the natural drainage pattern in a watershed.  A simplified example is show here (yellow area is HCZ, which lies partially in basin- the pink area):

I am having a lot of trouble coming up with a good method to map areas upslope and downslope of roads. My current methodology is based on @nagib's post:

Upslope and downslope distance or area using ArcGIS Desktop

This method may have promise, but details in those steps are sparse, especially for the all important and complex second step- path distance allocation.   When I tried to mimic nagib's method, I didn't get a very intuitive result. Not surprising, because I am not very familiar with calculating upslope/downslope using path distance allocation.  The result looked like:

In mapping HCZs it seems to me that accurately defining upslope and downslope areas of roads is a critical first step.  Any resources or thoughts on a very accurate methodology for getting these areas right?  Using ArcMap 10.4 primarily, but have access to QGIS and Whitebox GAT.

EDIT:

Per @nagib's request here is a third pic showing second pic with contours (white lines):

Also, here is a pic combining upslope/downslope areas, contours, pour point(orange point) at intersection of road(black line) and river(blue line), and the resulting subshed (white area) from that pour point:


texts: OpenLayers Geolocation not using View projection?

texts:
I'm struggling with Openlayer's Geolocation API, using OL v. 4.4.2. For some reason, geolocation.getPosition() returns the OL coordinates array in WGS84, even though the map's view projection is something else. Here are the relevant code bits:

const extent = [548576.0,6291456.0,1548576.0,8388608.0];
const center = [327305.0,6822495.0];
const proj = new ol.proj.Projection({ code: 'EPSG:3067', units: 'm', extent: extent });
const view = new ol.View({ center: center, projection: proj, zoom: 8 });

let geolocation = new ol.Geolocation({
  projection: view.getProjection()
});

let positionFeature = new ol.Feature();
let whereat = new ol.layer.Vector({
  source: new ol.source.Vector({
    features: positionFeature
  })
});

geolocation.on('change', function() {
  let coordinates = geolocation.getPosition();
  positionFeature.setGeometry(
   coordinates ? new ol.geom.Point(coordinates) : null);
});

let map = new ol.Map({
  layers: [basemap1, basemap2],
  target: 'map',
  overlays: [overlay],
  controls: [new ol.control.Zoom({className: 'custom-zoom'})],
  view: view
});

$('#track').change(function() {
  geolocation.setTracking(this.checked);
  this.checked ? map.addLayer(whereat) : map.removeLayer(whereat);
});

Looking at the browser console, geolocation tracking seems to work and also changes in the position are registered, but the location does not appear on the map, and I suspect this is due to the geolocation being still obviously in WGS84, unlike the map view and other layers.

Logging the geolocation, I get for example: [23.7633003,61.5019864], so obviously in WGS84, even though EPSG:3067 is intended. Any solutions or workarounds? I've thought of setting the geolocation projection initially to WGS84 and later on using ol.proj.transform to transform it to EPSG:3067, but I find this a dirty trick, since I figure the Geolocation API is supposed to do this with the "projection" property.


texts: Finding intermediate points and speed over a segment

texts:
I have a dataset which is in the form of ( Time , Car ID, Speed at that moment, Lat,Long) and I need to do some machine learning on the data to answer the following questions :

Best route from point A to point B
The speed at some point ( which can be on any segment) .

I am having a few problems because I have never dealt with this kind of data before ( I mostly worked with text before this ) :

The data is Hyderabad specific. I want to divide Hyderabad into segments but how do I go about dividing it based on points ? If I use 500m as minimum length and keep dividing ,  there can be many twists and turns which will go missing :(
If I get a random point , how do I decide which segment it falls on ? For example, I have two segments , each having start point ( lat,long) and end point ( lat',long') . If i get a third point  ,  (l1,l2), how do I decide which segment it falls under ?

I actually have more questions, but if these get solved first, it would be awesome.

Thanks a lot


texts:
I have a project on QGIS 2.18 which includes several raster and vectorial layers and need to export the full project as a high resolution Geotiff.

QGIS allows to export project as an image on print composer but it exports a low resolution image.

I am also able to export layer by layer but I need to export the full project as a single file.

Is there any plugin or approach that allows me to achieve that?


texts: QGIS upload Ortho NAIP raster layer problem - loaded a raster file, georeferenced, however is blurry

texts: Relatively new QGIS user. I am trying to correctly upload an ortho NAIP file I got from the USDA geospatial gateway. I successfully uploaded the raster layer from a .sid file of a county. As I zoom in to my area of interest, the resolution is quite blurry. What I need to do is draw plot polygons over top of the plot area in the center of this image, but it would be helpful if I can actually see the tractor tire tracks in between plots. Did I do anything wrong that created this blurry image, or is this as good as it gets?


texts: Raster calculator - problem with number variables in ModelBuilder

texts:
When I 'm calculating this expression:

("%Soil organic carbon (t/ha)%" / (10000 * "%Soil depth (m)%" * (Float("%Bulk density (kg/m3)%") / 1000))) * 100

Where:

"%Soil organic carbon (t/ha)%" and "%Bulk density (kg/m3)%" are raster layers

and

"%Soil depth (m)%" is any value or calculator expression with input values: 0.05, 0.15, 0.30, etc.

I 'm geting this error:

  ERROR 000539: Error running expression: rcexec()  Traceback (most
  recent call last): File "", line 1, in   File
  "", line 5, in rcexec   File "c:\program
  files\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 4483, in
  Times
      in_raster_or_constant2)   File "c:\program files\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper
      result = wrapper(*args, **kwargs)   File "c:\program files\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 4480, in
  Wrapper
      ["Times", in_raster_or_constant1, in_raster_or_constant2]) RuntimeError: ERROR 000732: Input Raster: Dataset
  0.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.050.05   does not exist or is not supported

  Failed to execute (Raster Calculator).

Оn the contrary, when I write down "Soil depth" by hand like this, everything is ok

How to fix this problem?


texts: Setting an item position (variable) from LowerRight?

texts:
I'm writing a standalone app with PyQGIS to export a pdf representing a building site. I managed to write the whole code to do it, but when I'm building my pdf, I want to set the position of my QgsComposerLegend object to the LowerRight angle of the QgsComposerMap object. The position of the QgsComposerLegend item has to be set progammatically because the layer set, and in fact its size, is variable (the program will be used for a couple of qgis projects...), so I need to set this postition using the lower right corner of the composerlegend object.

I found a way to do it:

    legend = QgsComposerLegend(c)
    legend.setTitle("Legende:")
    legend.model().setLayerSet([layer.id() for layer in layers])
    legend.updateLegend()
    legendSize = legend.paintAndDetermineSize(QPainter())
    legend.setItemPosition(x+w-legendSize.width(), y+h-legendSize.height(), legendSize.width(), legendSize.height(), QgsComposerItem.UpperLeft, False, -1)
    legend.setBackgroundColor(DEFAULT_BACKGROUND_COLOR)
    legend.setBackgroundEnabled(1)
    legend.setFrameOutlineWidth(0.5)
    legend.setFrameOutlineColor(DEFAULT_COLOR)
    legend.setFrameEnabled(1)
    c.addComposerLegend(legend)

In this code: layers is the list of the layers in the composermap, x y w and h are the coordinates (in mm) of the composermap for the pdf.

Here is the problem I get: when running my code, I have these messages and I don't want them to appear:

I would like to try something like this:

    # On met en place la legende:
    legend = QgsComposerLegend(c)
    legend.setTitle("Legende:")
    legend.model().setLayerSet([layer.id() for layer in layers])
    legend.updateLegend()
    legend.setItemPosition(newX, newY, newW, newH, QgsComposerItem.LowerRight, False, -1)
    legend.setBackgroundColor(DEFAULT_BACKGROUND_COLOR)
    legend.setBackgroundEnabled(1)
    legend.setFrameOutlineWidth(0.5)
    legend.setFrameOutlineColor(DEFAULT_COLOR)
    legend.setFrameEnabled(1)
    c.addComposerLegend(legend)

But I can't manage to make it work...
Could anybody help me ?


texts: Red marker cross on canvas can't be hidden or deleted

texts:
I've somehow created an indelible red X on my map in QGIS. I can't select it, and switching off all the layers doesn't hide it. It doesn't show in the print layout. Its a constant screen size, eg mm.
Is it some feature I accidentally switched on?
How do I get rid of it?
Running 3.14.6 on W10


texts: Publishing metadata (.xml) for every single row in table on GeoServer?

texts:
How do I publish metadata in xml format for every single row in a table via Geoserver?

In details, I know i can use GeoNetwork or OpenGeo Suite(Boundless) for publish metadatas but I am very curious about publishing the metadatas using just Geoserver. If I want to publish metadata files for single layers, I add a metadata links in description of layer in Geoserver.
Metadata Import in GeoServer for layer but it doesn't work for every single rows in table.

Note 1: I have already installed CSW plug-in for GeoServer.

Note 2: I cant find any examples metadatas which published in GeoServer.


texts:
I am trying to send (x,y) coordinates from script HTML code to PHP but it doesn't work:

-Script code:

searchControl.on("results", function (data) {
            results.clearLayers();
            for (var i = data.results.length - 1; i >= 0; i--) {
                results.addLayer(L.marker(data.results[i].latlng, { draggable: true, title: 'Local' }));
                x = data.results[i].latlng.lng;
                y = data.results[i].latlng.lat;

-PHP piece of code to get coordinates:

    <?php
$dbconn = pg_connect("host=*** dbname=*** user=*** password=***")
    or die('Not connected: ' . pg_last_error());

// Get Parameters

$x = $_GET['x'];
$y = $_GET['y'];

Is $.ajax the best option? How could it be possible with this?.


texts: Is it possible to make a geodatabase feature class invisible?

texts:
I am working on creating a report for a customer that requires using a couple tools in arcmap then export values to a dbf. They then have to run a crystal report that pulls data from Oracle in the cloud that we have granted restricted access.

I want to create a report directly in arcmap that would be incorporated into a tool, but I can't figure out how to access that data needed for the report without incorporating the restricted data (eg, using a fc with the data).  Is there a way to have a feature class embedded in a gdb but restrict the access to that fc so that the customer doesn't see it?

I understand this is not ideal, but it is kind of a bandaid till we get the funding to build out our online tool.


texts: Check if feature attribute is empty in QGIS with Python

texts:
I want my Python script to do something when two feature attribute values are None at the same time.

So I made this if statement: if feature["Name"] and feature["Zon"] is None:, but the code stops there and doesn't continue. Before the if statement I print both feature attribute values and get the following: (None, None). So the code should theoretically work, but it doesn't. What am I missing?


texts: Is it possible to merge different vector types in qgis

texts: I would like to merge a point vector layer with a polygon layer. There are some attributes in the polygon layer which I would like to be in the point layer!


texts:
I've been asked by my manager to come up with some sort of competency measure for an employees understanding of scripting/programming. Because I am a self taught programmer, I really don't have a formal understanding of how to measure ones level of understanding when it comes to script writing.

We have been using Python in ArcGIS and QGIS. We use VB.Net with our MS Office products and we have started to dabble in Java for web mapping.

I'm wondering if any of you GIS managers out there have some sort of matrix/method to determine your employees skill level when it comes to scripting.


texts: Addressing based on centrelines, or assigned at parcel creation?

texts:
I'm still pretty green when it comes to GIS and using ArcGIS products, so hopefully this question doesn't sound too dumb.

I am just wondering what the best way to go about assigning addresses to newly created parcels within a new subdivision.

Should I create the road centrelines and add in the FromLeft, ToLeft, etc. attributes, and assign addresses based on this, or should I assign an address when the subdivision is created, and update the centreline attributes to correspond with the created parcels.


texts: Disable selection of clusters OpenLayers 3

texts:
At the moment I have my OpenLayers 3 map set up to cluster icons upon zooming out. I want to disable the ability to select clusters but allow the user to select individual markers if they zoom in. Preferably, clicking on a cluster could select each individual feature in the cluster (in their correct locations). My select code is like so:

var selectClick = new ol.interaction.Select({
    condition: ol.events.condition.click,
    style: selectedMarkerStyle,
    addCondition: function(mapBrowserEvent){
        var features = null;
        map.forEachFeatureAtPixel(mapBrowserEvent.pixel, function (feature_, layer) {
            features = feature_.values_.features;
        });
        if(features !== null && features.length > 1){
            return false;
        }
        return true;
    }
});
map.addInteraction(selectClick);

However the addCondition doesn't seem to be working. At the moment clicking on a cluster selects one of the features in the cluster, in the location of the cluster, rather than the features correct location.

Is there any way to fix this?


texts: Converting GeoJSON to ISOXML ISOBUS ISO 11783-10

texts:
How do I convert some GeoJSON files to ISOXML files (ISO 11783-10)?

I found several questions on this topic:

Converting GeoJson/Shp to Iso-xml (ISOBUS ISO 11783-10)
https://stackoverflow.com/q/48623806/13070285

I am trying to contact the author of these questions directly to find out more details.

But maybe someone has experience converting files between these formats.

Especially using the ISOXML plugin for ADAPT framework:

https://github.com/ADAPT/ISOv4Plugin
https://github.com/ADAPT/ADAPT

I try to import an ISOXML file into an ADAPT Data Model using the ISOXML plugin, understand the data structure inside the model and try to export back to the ISOXML format.

https://github.com/ADAPT/ISOv4Plugin/blob/bc5ff975338abe85963f882945a7ed61e7ac470d/ISOv4Plugin/Plugin.cs#L28-L42

I am stuck on the fact that I can’t understand if I can export the model back to ISOXML file. I asked these questions directly to the developers on GitHub and by email, but haven't received an answer yet.

Maybe someone was dealing with the ISOXML format or tried to do the same with the ISOXML plugin for the ADAPT framework.


texts:
I have six polygons, all of (fairly) similar size. The first one is 8 acres, the others are obviously wrong with 120,000 acres, and 106,000 acres and so on. I am calculating with the field calculator, Geometry:

area* .000247

Is there something I should check? Is there another method I could use?


texts: Using ArcMap, I am trying to classify/consolidate the symbology of these agriculture classes into 3 categories: Hay, Pasture/Grass, and Cropland. I also need to find the sum of acres of each of these 3 categories after reclassifying them. The file is currently in a raster format.


texts: Automatically dividing multiple polygons into halfs

texts:
I have this shapefile with multiple polygons (equally sized hexagons). My wish is to divide all hexagons by two. The goal is then to have two layer, one layer with all the left half of the hexagons and the second layer with the other half. If it is left, right, up or down doesn't matter really.

The difficult part for me is to find an automated mechanism that can handle this.

The purpose for this procedure is that I want to display two different values for the same hexagon. Therefore, my idea was to separate it like this.


texts: FeatureCollection feature.byFeature histogram or column chart bucket size

texts:
I am working with a featureCollection of river data which is segmented by river reaches (there are hundreds within a few km). This collection includes data of reach length and reach discharge. I am trying to plot a histogram or column chart where the x-axis represents reach length and the y-axis the reach discharge.
So far I've managed to build a column chart but where length is represented by individual columns rather than grouped by buckets,

or to build a histogram for reach length and reach discharge separately.

What I would really like is to be able to group hAxis values into buckets for the first chart, but I don't know if it is possible with GEE:

This is the code I've written (studyArea can be any random polygon):
var dataRivers = ee.FeatureCollection("WWF/HydroSHEDS/v1/FreeFlowingRivers");
var dataRiversStudy = dataRivers.filterMetadata('INC', 'equals', 1).filterBounds(studyArea);

print('dataRiversStudy:',dataRiversStudy);

Map.addLayer(dataRiversStudy,{color: "darkblue",width: 2.0},'Rivers');

var chart = ui.Chart.feature.byFeature({features:dataRiversStudy, xProperty:'LENGTH_KM', yProperties:['DIS_AV_CMS']})
                .setChartType('ColumnChart')
                // .setOptions({bucketSize:2,maxBuckets:10}); // This does nothing

print(chart);

var chart2 = ui.Chart.feature.byFeature({features:dataRiversStudy, xProperty:'LENGTH_KM', yProperties:['DIS_AV_CMS']})
                .setChartType('Histogram');

print(chart2);

And here is the link to the code: https://code.earthengine.google.com/4cc589ee05c037494eb20839402096da
Thanks!


texts: Need to resize XY data point symbols based on numbers of sightings in a third column of data. Was fairly staight forward in ArcGIS. Can't figure it at all in QGIS 3.0. Most answers I've seen refer to a graduation tab in the Style Menu in Properties but I can't even locate the style menu shown for earlier versions of the program let alone a Graduation Tab.


texts:
I am trying to implement the code found at this site: http://geospatialpython.com/2015/05/clipping-shapefile-in-pure-python.html
The function I am trying to implement:
import shapefile

# Create a reader instance for our US Roads shapefile
r = shapefile.Reader("roadtrl020")

# Create a writer instance copying the reader's shapefile type
w = shapefile.Writer(r.shapeType)

# Copy the database fields to the writer
w.fields = list(r.fields)

# Our selection box that contains Puerto Rico
xmin = -67.5
xmax = -65.0
ymin = 17.8
ymax = 18.6

# Iterate through the shapes and attributes at the same time
for road in r.iterShapeRecords():
    # Shape geometry
    geom = road.shape
    # Database attributes
    rec = road.record
    # Get the bounding box of the shape (a single road)
    sxmin, symin, sxmax, symax = geom.bbox
    # Compare it to our Puerto Rico bounding box.
    # go to the next road as soon as a coordinate is outside the box
    if sxmin <  xmin: continue
    elif sxmax > xmax: continue
    elif symin < ymin: continue
    elif symax > ymax: continue
    # Road is inside our selection box.
    # Add it to the new shapefile
    w._shapes.append(geom)
    w.records.append(rec)

# Save the new shapefile! (.shp, .shx, .dbf)
w.save("Puerto_Rico_Roads")

I am getting the following error: AttributeError: 'Writer' object has no attribute '_shapes'
I think this issue may be connected to the issue discussed here: https://github.com/GeospatialPython/pyshp/issues/113
However when I convert the lines:
w._shapes.append(geom)
w.records.append(rec)

to something like:
w.shape(geom)
w.record(rec)

This does not help. And I am unsure of how to append multiple records or shapes to the writer.


texts: I would like to add a google base layer to my openlayers map but it load correctly

texts:
This is my script at the top of my HTML file

<script src="http://www.openlayers.org/api/OpenLayers.js"></script>
<script src="http://maps.google.com/maps/api/js?v=3&amp;sensor=false">
</script>
<script type="text/JavaScript" src="map.js"></script>

This is my .js

var map;
function initMap() {
//initialise variables
 var osmlayer;
 var trackloca2;
 var bmenatrack2;
 var pits3;
 var boundary3;
 var track2;
 var pits2;
 var boundaryline;
 var gmap;
//create new map object
map = new OpenLayers.Map("mymap");
 //create new OSM base layer
 osmlayer = new OpenLayers.Layer.OSM("Open Street Map");
 gmap = new OpenLayers.Layer.Google(
    "Google Streets",
    {numZoomLevels: 20}
    );

I have added gmap to the addlayers at the bottom but when I switch from the open street map to the google layer all I am getting is a white screen with my geoserver layers showing,


texts: Create fishnet of irregular polygon study area

texts:
I would like to create a fishnet (grid of 150x150m cells) for a study area (irregular shaped polygon). However, the Create Fishnet tool creates a fishnet representing the full extent of the study area (max/min long/lat) and I only need grid cells within the study area.

I can't seem to select by location (freezes) and then create layer from selected features that intersect the study area polygon. The # of grid cells is possibly too large (study area is ~440,000 km2, fishnet includes 30 million features).

Is there some other way to limit the Create Fishnet tool to only create polygon grid cells within the study area boundary?