texts: Getting distance between 2 points using GeoDjango?

texts:
I have 2 locations defined by gps coords, lat / long like returned by google maps:

http://gmaps-samples.googlecode.com/svn/trunk/geocoder/singlegeocode.html

I need to calculate the distance between them. I know I could use the google API but I'll be processing bulk queries so I'd rather do it on my own server.

I've spent a few hours with the docs, installed geodjango OK, but I can't seem to find an example of this. Everything in the literature is way more complicated than I need.


texts: Can I find out the exact areas which each category covers?

texts:
I have accessed a land use map and clipped the appropriate area I need to analyse. Is there a
way I can find the exact areas which each category covers?


texts:
I'm a new user in Leaflet and I create a weather map with a lot of weathers symbols around the world.
I can already filter the number of symbols on the map with a function :

var filteredSites = filterSitesBySpacing(sites, 0.8);

But I would like to filter according to the level zoom. So, I use if condition :

map.on('zoomend', function() {
    if (map.getZoom() <10){
  var filteredSites = filterSitesBySpacing(sites, 0.8);
    }
    if (map.getZoom() >= 10){
  var filteredSites = filterSitesBySpacing(sites, 0.8);
    }
 }

But, it doesn't work...
Could you explain me what is the error in my "if" condition ?


texts: Create a rectangle with specific height from two points

texts:
I need to be able to create a rectangle with height x from two points. This is probably not clear enough so I tried to draw it:

The two points that I have are the centroids of those circles, and what I need is the rectangle in orange. I have little experience with GIS, and my idea is probably silly, but I'm thinking that if I create a line connecting the two points, and find the intersection in the circle and then find a point 90 degrees from that point, I could calculate the 4 points that I need to create a rectangle.
This example in real life is like a tractor working from point A to point B, and there's this equipment attached to the tractor which has a known width, and the rectangle tells me the area that this tractor did something.
I'm working with this 'map' on a web page, in JavaScript with Turf.js.


texts: Use full QGIS functionality on Server with GUI

texts:
I would like to know if there is an easy and efficient way to use the full QGIS desktop experience on a server.

Currently I installed QGIS desktop on an Ubuntu Server and I am able to open the graphical interface over SSH:

ssh -X username@remoteserver

and then simply opening QGIS by typing "qgis" which runs the binary on the server but opens the gui on my desktop pc (Ubuntu as well).

It has full functionality but the interface is not very efficient (long rendering) if I access the local network from outside (from inside it is fine though).

I saw information on QGIS server but it looks like a stripped down WMS version of the original QGIS. I wondered if there was something that offers me a broader range of function with the efficiency of QGIS server technology for remote usage...?

Why do I need full functionality? Well, instead of only displaying and filtering data (think of WMS) I would like to use this setup to create spatial models that require QGIS algorithms. Those models can often take up long processing hours and by executing them on a server I could free my Desktop PC from the processing overload.


texts: Esri Geoportal Server on Linux install error?

texts:
I'm following the How to Set Up Esri Geoportal Server on Linux white paper.
Going everything OK until I have to run ./grants_linuxpg.sh (page 10)

./grants_linuxpg.sh: 20: ./grants_linuxpg.sh: Syntax error: "(" unexpected

Did someone already had this error? What do I wrong? I just can't believe that there is a syntax failure in the original script downloaded from the geoportal sourceforge website?!


texts: How can I have the centroid of each chunk of a catalog?

texts: I use lidR package. I have a catalog and I want to save as a shapefile the centroid of each chunk make with opt_chunk_size of a catalog ?


texts: How to change sign of z value in a large dataset?

texts: I have data in an XYZ format (saved as a .txt file) that I am importing into ArcGIS and then converting to a raster using the point to raster conversion. However the Z field in the data represents a depth below the sea surface, but has been output as a positive value by the source software. In order to plot this correctly in ArcGIS is there a simple way to convert these values to negative values or to plot them as subsea? The table contains almost 40,000,000 rows so I cannot do the conversion in a text editor. I have tried using the Field Calculator function in the table, however this option is grayed out?


texts: Where can I find GIS data to test GIS coordinate operations?

texts:
I am trying to find some data that can be used to test GIS coordinate operations from one CRS to another CRS. So far all I have found is the Gold Data set and also various example points in the EPSG G7-2. The gold data is OK for testing some projections from WGE and the few points in the G7-2 document are good for testing individual operations. What I need now though are some points that can test implicit as well as explicit concatenated coordidnate operations such as EPSG:4277 to EPSG:4230 for example, which I had found a sample for in http://www.epsg.org/Exchange/1065.pdf. There has to be a good source out there with more than one point and CRSs using multiple datums. Fly, Google monkeys, fly!

EDIT:

I also found something called GIGS


texts: What coordinate system is this netCDF raster?

texts:
I obtained a netCDF raster file, but I couldn't get any metadata to get the name of the coordinate system the raster has been built on. The raster itself doesn't have any coordinate system embedded. I thought it would just be a WGS84, and it looked like that at first glimpse, but with further investigation in ArcMap, I saw that it is a rather not common system. Here is how it displays:

The orange raster is a normal raster in WGS84 which I have inserted here for comparison purposes. The purple one is the raster with the unknown coordinate system. Do you have any clue what this might be?

Some updates: Here is the netCDF raster: https://www.dropbox.com/s/nottbl9yt6dwss6/sic_average_nclimate.nc?dl=0
I was also able to get some metadata from the image provider:

netcdf sic_average_nclimate {
dimensions:
    nlon = 361 ;
    nlat = 90 ;
    nseas = 4 ;
variables:
    float SIC_Change(nlat, nlon) ;
        SIC_Change:Title = "Gridded Multi-Model Ensemble Mean Annual Mean Change in Ice Concentration 21C-20C" ;
    float SIC_Season_Change(nseas, nlat, nlon) ;
        SIC_Season_Change:Title = "Gridded Multi-Model Ensemble Mean Seasonal Mean Change in Ice Concentration 21C-20C" ;
    float SIC_Change_STD(nlat, nlon) ;
        SIC_Change_STD:Title = "Gridded Multi-Model Standard Deviation of the Annual Mean Change in Ice Concentration 21C-20C" ;
    float SIC_Season_Change_STD(nseas, nlat, nlon) ;
        SIC_Season_Change_STD:Title = "Gridded Multi-Model Standard Deviation of the Seasonal Mean Change in Ice Concentration 21C-20C" ;
    float LAT(nlat) ;
        LAT:Title = "Latitude" ;
    float LON(nlon) ;
        LON:Title = "Longitude" ;

// global attributes:
        :Title = "Ice Concentration metrics for Model subset as in Figure 1 of NCLIMATE paper" ;

They show the boundary lat-long, but apparently not any information regarding the coordinate system.


texts:
The computer program that I'm working with requires an 8 digit UTM Easting for the coordinates of shapefiles (i.e., the UTM zone number must be included with each coordinate within a shapefile).  The program was designed in Europe and considers all of my shapefiles to be extremely far away, due to the two "missing" digits.

I'm using ArcMap 10.1 and also downloaded Quantum GIS, as the individual I'm working with overseas suggested that QGIS could add the UTM zone number to the coordinates.  I haven't found any way to do so, so if anyone has any ideas, I'd be very grateful.

Thanks in advance!


texts: GRASS v.rast.stats missing in SEXTANTE toolbox

texts:
I am new in the QGIS world and also with GRASS. I am really interested in using v.rast.stats as it is doing exactly what I need. My problem is that I am not able to find it inside SEXTANTE in QGIS and I don't now how to make it appear.

Can anybody help me?


texts: How can I extract time and elevation from this data?

texts:
I download metadata from
http://geoportale.lamma.rete.toscana.it/MapStore/public/ based on this tutorial
http://geoserver.geo-solutions.it/multidim/en/accessing_multidim/rtx/wind_arrows.html and extract it to a folder, so there is some .tiff files in it.
then I opened my Geoserver 2.8.1 and on this path: Stores->add new Store-> ImageMosaic ( image mosaicking plugin) then choose the folder contains tiff files.
then in the layer configuration, dimension tab, I should see that the time and elevation attributes are enabled as it is mentioned in the tutorial but they are not and The following message appears next to each of the respective enable checkboxes:

Time

Cannot enable, no attribute of type Date found

Elevation

Cannot enable, no attribute of type Number found

I already saw this question How to add time and elevation dimension attributes to Geoserver rasters? but that is not my case since I have not insert Geotiff Data into my Postgres or whatsoever.
I spend two days on solving this and read a lot about it with googling the issue.


texts: Check if a VERTEX of one polygon falls anywhere within the boundary of another polygon (in Python)

texts:
This is different than simply containing, overlapping, etc. For example, imagine a circle partially overlapping a square:

The circle and square both overlap each other, both intersect each other, but only the circle has vertices within the square. (The square only has four vertices, none of which are inside the circle. In practice, since I'm dealing with convex hulls that are creating arcs, the "circle" would be a high-N N-gon).

How can I check for this condition? I was hoping one of Shapely's polygon methods would help, but they don't seem to answer this question. I also tried using difference in both directions and seeing which one contained all the points of the other, but this not only is very slow but is prone to making accidental holes in the output shapefiles later.

I plan to ultimately end up with a result like this, in case it is helpful:

but it'll be done to a dynamic number of dynamically placed centroid-esque shapes


texts: QGIS - Tree crown polygon maping with Geometry generator

texts:
I am new to QGIS and I have stumbled on this automatic way to map tree crowns using 4 points - north, east, south, west.
Most informative post was this and Solution 2 (8-vertex polygon).
The ellipse method works like a charm, but the drawn area is to vague for my purpose.
There I found code for exactly what I need:
make_polygon(make_line(
 translate($geometry,0,max("N",0.2)),
 translate($geometry,0.6*max("E",0.2),0.6*max("N",0.2)),
 translate($geometry,max("E",0.2),0),
 translate($geometry,0.6*max("E",0.2),-0.6*max("S",0.2)),
 translate($geometry,0,-max("S",0.2)),
 translate($geometry,-0.6*max("W",0.2),-0.6*max("S",0.2)),
 translate($geometry,-max("W",0.2),0),
 translate($geometry,-0.6*max("W",0.2),0.6*max("N",0.2)),
 translate($geometry,0,max("N",0.2))))

Problem is that no matter what I try, output is still NULL. Even if I change NESW to numbers, its sill NULL. It should work, at least it has worked in the past for the person who wrote the code.
Naming for NESW is correct. It was the first thing I checked.
I have a general idea how it should draw the polygon.
That post was made few years ago. Has something changed?


texts: PROJ4 / PostGIS transformations between WGS84 and NAD83 transformations in Alaska

texts:
I'm trying to implement the recommended WGS84 - NAD83 transformation for Alaska using PostGIS.  This post: Geographic Transformations in Alaska Between NAD83 and WGS84 indicates the proper transformation in ArcGIS to be WGS_1984_(ITRF00)_To_NAD_1983.  I haven't been successful at finding this in spatial_ref_sys.  I tried adding my own transformation into spatial_ref_sys using the following proj4text:

+proj=longlat +ellps=GRS80 +datum=NAD83
+towgs84=-0.9956,1.9013,0.5215,-0.025915,-0.009246,-0.011599,-0.00062
+no_defs

The towgs84 values were copied from the listed settings that ArcMap shows when you choose this transformation (reversing the sign because it's to WGS84 in PROJ4 and from WGS84 in ArcMap).  In PostGIS, I use it like this, transforming WGS84 to NAD83 using my custom projection, then transforming it again into 4269 (or whatever NAD83 projection I really want, like 3338).

ST_Transform(ST_Transform(geom_wgs84, 900999), 4269)

The problem is that the transformed values don't seem to match the transformations that ArcGIS are using, so if I have a mix of WGS84 data and custom-transformed into NAD83 data, ArcGIS's WGS84 -> NAD83 transformation is different from PostGIS and the data doesn't line up.  The errors are small (1-2 meters), but large enough to see.

Can anyone suggest an appropriate PROJ4 string to replicate the results of the ArcMap transformation, or something better for the state of Alaska such as using a NAD grid?  Or tell me what I'm doing wrong with my custom projection?


texts: Subtracting point values from raster maps in QGIS

texts:
I have point values (x,y,z) and a raster map with z-values. I would like to subtract the z-values of the raster map from the z-values of the points.
The resulting data has to be point data (x of point values, y of point values, z of point values minus z of raster map values).
How do I do it? By hand won't be an option because there are way too many points. I use QGIS.


texts:
I want to select many objects on the screen and select other objects that share the same attributes with the chosen group.

So I'm creating a tool to select many object based in five tables' columns attributes.
The tools works fine but when I select 3 or more objects the query became much bigger and I lost part of them.

Is there a way to do a query using to variables?
Something like: "Run command varA + varB"
Or is there a way to increase the string's size?


texts: I'm wondering if anyone has heard of or come across any examples of scanned and georeferenced Sanborn Fire Insurance maps being served as a WMS.  I know there are copyright issues, so presumably only for maps before 1922.  I'm curious about the possibility of doing this myself...  Thanks!


texts: Clipped raster excludes some pixels at the edge-QGIS

texts: I am trying to clip a raster using mask vector polygon layer. But clipping it causes some pixels at the edge to not get included in the output raster. I want to get all the pixels in and the ones which fall on the line to be included into the outptut file raster. I am using gdalwarp to clip the file something like: gdalwarp -q -cutline Masklayer.shp -tr 0.125 0.125 "input.tif" "output.tif"


texts: Sharing QGIS Project that connects to local postGIS database

texts: I am looking to share a QGIS project with a few others outside my network. The QGIS project is using layers from a postGIS database I created locally. I am wondering what the best way is to go about sharing the project and granting access to the database files. Is it possible to allow them remote access to my database, or should they just create their own databases with the same files?


texts:
The images below is a snippet of 2 shapefiles.

1) every building in the UK (13,971,598 features)

2) blue polygons

The goal is to clip or select all the buildings that fall within the blue polygon features. I'm wondering what method would be most time efficient. I have tried select by location but ran into the error posted below. Other than using GDAL/OGR or SAGA algorithms to clip the buildings layer which I imagine would take >24hrs, would there be an alternative method? I would be interested to know how you would tackle this so that I may test it.
My experience of working outside of QGIS desktop is very limited so please go easy on any jargon that you might find someone like myself to be unfamiliar with.


texts: Can any WFS server accept ECQL filter as a POST request parameter

texts:
I believe GeoServer requires that filters passed in a WFS POST request are in the OGC XML based format.

If you make a GetFeature request to GeoServer using a GET request the filter can be CQL or ECQL and GeoServer will accept it.

Will GeoServer ever be able to accept CQL in a POST filter parameter or provide an endpoint which will convert it on the fly into the XML based filter?

Alternatively, do any of the other web servers handle this better?

Currently our layers are filtered using CQL for the WMS requests so it would be great if the WFS GetFeature requests were using the same filter format.


texts: osm2pgsql style - removing amenity=school but keeping the rest of amenities

texts:
I am trying to import a PBF file using osm2pgsql, and the file looks like this:

...
node,way   amenity             text         polygon
# remove the school amenity data
node,way   amenity=school      text         delete
...

The problem is the amenity=school tags still make it in the database, is there a way to include all amenities but exclude the amenity=school tag, or should I just process it via the SQL DELETE statement in postgres?

I understand the delete flag prevents the column from being created, so this might be the reason. What are the possibilities in this use case?


texts: How to set Scale Factor for Lambert Conformal Conic (LCC) projection in ArcGIS Desktop?

texts:
The projection for my state according to my country's Mapping Agency is in a specific Lambert Conformal Conic (LCC) projection. They have given values for the following parameters:

False Easting
False Northing
Central Meridian
Standard Parallel 1
Standard Parallel 2
Scale factor
Latitude of origin
Datum/GCS

When I go to ArcCatlog, to update an existing LCC projection, I can see all the parameters except for the Scale Factor:

How do I enter the Scale Factor? The ArcGIS help mentions this parameter but does not mention how it to be used.

A semi-related question is: How important/required is the Scale Factor? I have seen several LCC projection without the Scale Factor. Even Snyder's Map Projections: A Working Manual does not mention it, when describing the Lambert Conformal Conic projection.


texts: Removing same name layers from layer table in geotool?

texts:
I am making a geotool application in which I read the data from postgresql database. In this application if i make a buffer on a selected polygon then a buffer is created and buffer layer shown in layer table of JMapFrame. After that i make another buffer on another polygon then there are two layers with same name  shown in the layer table in JMapFrame and both the layers shows the buffer on last selected polygon.
How can I replace the same name layer with the new one from the layer table in JMapFrame.

My code runs on the click of "Buffer" Button

And Result on Jmapfrmae

In this JMapFrame all the khasra_buffer layers are the same i.e. the current buffer which is created at Id 80, Thus replacing previously created two buffers.

I want that old khasra_buffer layer should be removed and the only buffer layer that should be shown is of the current buffer on JMapframe.


texts:
I'm trying to use Animate Columns in MMQGIS plugin. I don't have a lot of input as of how to use it, apart from the instructions here: http://michaelminn.com/linux/mmqgis/

What I don't understand is how to calculate the offset values. The instructions say:

  Offsets for points may be calculated by:  offset_long = (end_long -
  start_long)   offset_lat = (end_lat - start_lat)

Does that mean I should have a column subtracting the starting Lat from the ending Lat etc.? eg. 2538295.09836-2544008.68716 = -5713.5888099. But that doesn't make sense to me. Alternatively, I tried having a text column with the 2538295.09836-2544008.68716. But the tool, doesn't recognize text columns.


texts: Is there a location where one can download all of the gpx files for OSM?

texts: Where can I download all the gpx file in OSM.org website? Is there a package include all the gpx file? I don't want to download these gpx files one by one...


texts: SAR - Using the Intensity or the Amplitude for the classification

texts:
In a Sentinel-1 Level 1 Ground Range Detected SAR data product, Intensity = Amplitude * Amplitude.

Should the intensity image or the amplitude image be used for radiometric calibration when using the ESA Sentinel-1 Toolbox software? The ultimate purpose is the classification of geophysical features, e.g. land, water, etc.


texts:
I'm trying to install the i.atcorr extension on GRASS 7.4.0 using the following (Running Linux Mint 18.3):

g.extension extension=i.atcorr url=trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr

Which always returns

ERROR: Extension <i.atcorr> not found

The following also returns the same error message:

g.extension extension=i.atcorr url=http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr?format=zip

Is this the correct way to install extensions from OSgeo, or is this even the correct url to the extension source code?


texts: Is it possible to delete ArcGIS imageserver using a script?

texts:
I have an app that allows users to upload TIFF files so that an ImageServer is created in ArcGIS 10.5 Server. However, I don't want to keep the images over a day.

Is there a way of automatically removing an image service using ArcPy?

The code below does not work as it seems to be designed for Desktop:

arcpy.Delete_management(url)

The same question is asked in Esri community here and here with no answer.


texts: Dynamic date range for layout element in ArcMap 10.5

texts:
In ArcMap I would like to setup an MXD to have a dynamic element that displays a date range.  It is easy enough to add dynamic text that displays the current date using the following:

<dyn type="date" format="short"/>

However, I really need to have a dynamic element that displays the a date range for the last 8 months.

In other words, if today is 10/23/2017 is it possible to create a dynamic text tag that would display "2/23/2017 through 10/23/2017"?


texts: River channel is cut after rendering flow accumulation in hec-hms?

texts:
I'm trying to generate flow accumulation using hec-geohms and I have already done some terrain processing using the DEM.
I was really confused what happened, why has the river channel been cut?

what could be the possible reasons why this occurs?
could this be a bigger error in my hydrologic model?

Here are the screen shots of what has been generated.
DEM overlain with river networks

Flow direction

Flow accumulation

Flow values


texts: Landscape change metrics for vector data

texts:
I am new to GIS.
I am currently working on a problem which involves two vector datasets about habitat and land use ("now" and "future"). I need to generate metrics like "number of patches", "mean patch size", "edge density", etc, to complete the analysis of landscape change. I see people using Fragstats a lot but I have vector rather than raster data.
Is there any tool or plug-in available to generate these from vector data?
I am using ArcGIS Pro 2.6.0.
I found one paper from four years ago talking about Arc_LIND, but I cannot find this plug-in anywhere.


texts: Local Copy of OpenStreetMap of my City in PostGIS Database

texts:
I would like to use OpenStreetMap vector data for land cover analysis in my city.
Here is my goal:

I set up a PostGIS database with the city area extent (around 20x20 km) on my server.
I will download the polygon data from OpenStreetMap. Maybe I could use a command-line tool like osm2pgsql
The local copy will be "read-only".
On a daily basis or more often, I want to "synchronize" my local PostGIS copy of the OpenStreetMap data of my city with the latest online edits to ensure that my data is up-to-date.

Is this possible without re-downloading all of the OpenStreetMap data from my region every time? (I would like to download the 'latest changes' only and somehow merge them with my local PostGIS OpenStreetMap copy)


texts: ArcGIS Select by Attribute SQL syntax problem

texts:
I have a list of aprox 150 unique road names in Excel. I copy and paste (as text option) that list into Word (2013 version). I then "replace" all the end paragraphs (^p) with quotes and commas (','). This gets the list to go from:

Cambie
Blanchard
Cook

to looking like

'Cambie','Blanchard','Cook'

Then within Arc's Select by Attribute I use this text in the query

Road_Name IN ('Cambie','Blanchard','Cook')

which gives me an SQL error. Any idea why I would be getting the error?

From trial and error I have come to realize that when I paste the text the quotes are a different length than if I had just typed them. See screen shot - the first road name's quotes (which I pasted) are shorter than the second road name's quotes (which I typed out in the query box).

When I query the one with the longer quotes by itself it works. Why would my pasted quotes be shorter? Any suggestions as to how I can query the 150 road names without having to re-type every quote?


texts: Read M value (CoordinateXYZM) from shapefile with geotools

texts:
After successfully writing M values into a shapefile (.shp, only points for now), I now also want to read the generated file and print the full coordinate set (x, y, z and m) using geotools (23-SNAPSHOT).

I found two different ways to do it (see code below) but both only give me a normal Coordinate object that has an M attribute, which is always NaN. I know that writing the file was successful because QGIS can see the full coordinates (including M).

How do I read the file, while still getting the actual M values?

Code I've tested:

1.From official website:

File f = new File(path);
Map<String, Object> map = new HashMap<String, Object>();

try {
    map.put("url", f.toURI().toURL());
    DataStore dataStore = DataStoreFinder.getDataStore(map);
    String typeName = dataStore.getTypeNames()[0];

    FeatureSource<SimpleFeatureType, SimpleFeature> source = dataStore.getFeatureSource(typeName);
    Filter filter = Filter.INCLUDE; // ECQL.toFilter("BBOX(THE_GEOM, 10,20,30,40)")
    FeatureCollection<SimpleFeatureType, SimpleFeature> collection = source.getFeatures(filter);
    FeatureIterator<SimpleFeature> features = collection.features();

    while(features.hasNext()) {
        SimpleFeature feature = features.next();
        List attributes = feature.getAttributes();
        Collection properties = feature.getProperties();

        if(attributes.get(0) instanceof Point) {
            Point p = (Point) attributes.get(0);
            Coordinate c = p.getCoordinate(); //"c.getM()" returns "NaN"
            CoordinateXYZM cxyzm = new CoordinateXYZM(c); //"cxyzm.getM()" returns "0.0"
        }
    }
} catch (Exception e) {
    System.out.println("Exception: "+e);
}

2.From a question here on Stackexchange:

File f = new File(path);

try {
    FileDataStore store = FileDataStoreFinder.getDataStore(f);
    SimpleFeatureSource featureSource = store.getFeatureSource();
    SimpleFeatureCollection featureCollection = featureSource.getFeatures();
    SimpleFeatureIterator it = featureCollection.features();

    while(it.hasNext()){
        SimpleFeature feat = it.next();
        Geometry g = (Geometry) feat.getDefaultGeometry(); //"org.locationtech.jts.geom.Geometry", NOT "org.opengis.geometry.Geometry"!
        Coordinate c = g.getCoordinate(); //"c.getM()" returns "NaN"
        CoordinateXYZM cxyzm = new CoordinateXYZM(c); //"cxyzm.getM()" returns "0.0"
    }
} catch (IOException e) {
    System.out.println("Exception: "+e);
}


texts: Advice on choosing ideal software for consultancy

texts:
I am taking the step to start a small business consultancy using digital map technology to assist retail business to better understand their environment and make informed location based decision. My focus is retail business and service planning, however, my challenge is the choice of software.

ESRI business analyst is not available to the part of the world where i live due to inadequacy of data but i don't want this to deter me from starting the business. I know there is always a way around obstacles if you work at it long enough.

While i work on solving the data adequacy issue i want to know if map-info professional will be a good software to purchase (I've never used it before). Any advice?


texts: Importing Functionality from Python Addin into script

texts: I am not sure if what I am trying to do is possible without using Arcobjects, but what I want to do is use a piece of a prebuilt addin created by ESRI and then add it into my own Arcpy script. The specific addin is the Curves and Lines addin and the Curvinator functionality is what I am after. I know it is possible to call a script from within another script but I can't figure out how to call the toolbar button from within a script or whether it is even possible. The addin is available at http://www.arcgis.com/home/item.html?id=b59665cb906a435da3f8c45332d7cfdf.


texts: Importing and converting coordinates to ArcMap and calculating distance between them?

texts:
I have a table in excel format with the following two columns:

Now, I need to import this table to ArcMap and calculate the distance between the two points for each row and write the calculated distance (kilometer) automatically in a new column "C".
Firstly, I use the tool “Excel to table”, but then:

What is the quickest way to convert the two columns in a format to use the “x,y, coordinate tool”? and has ArcMap a tool to calculate the distance between the two points of each row, when I do not integrate the whole street network?

I guess that I need to iterate through the table, but I have no straightforward idea which tools are the best to use.

I am working with ArcMap 10.5.1.


texts: Finding relation between Nominatim row and gb_postcode table row?

texts:
Currently I have PgAdmin 3 open with a copy of the Nominatim DB I processed on a Linux virtual machine. I imported the gb_postcode.sql database successfully from the Nominatim documentation but I'm struggling to find a way to match a row in my main placex table in Nominatim to my gb_postcode table that's been imported.

Both tables provide a Geometry row, of which I tried using ST_Intersection on but the result was beyond inaccurate (returned 1.6m records rather than a single one, or a small subset of data).

Is there any way currently that will allow me to find a relationship between a Geometry object in a Nominatim record (such as 0102000020E61000000500000014A5CEED6F5A01C07D6F78E68A1A4A4091E398767B5A01C08A517C21891A4A404221020EA15A01C02E05FFB6821A4A405698631E9E5A01C0CBC1C7BB7E1A4A40A53D14BB6C5A01C026A36F777B1A4A40) and the Geometry object in the gb_postcode table (such as 0101000020E61000005D96EA24715B01C0A76FB0828F1A4A40)?

For examples sake, these two Geometry values are matching for the same street in the UK, but I'm struggling to relate them programmatically.

Sample postcode data:

ID - Postcode - Geometry

627881;"HG4 3NJ";"0101000020E610000023F05E6B62D4F9BF9C759A1C19134B40"
627882;"HG4 3NL";"0101000020E6100000EAA10AA1EAEEF9BFBEA275E3EF134B40"
627883;"HG4 3NN";"0101000020E6100000AE7E4727E1F1F9BF4AE5506C06134B40"


texts: Getting complete address information from OSM without Geocoder

texts:
I know a question like this is asked very often. But the answers did never really satisfy me.

I want to extract addresses and maybe POIs from OSM. Nothing more. I am not interested any additional information.

So I inspected some OSM raw data and found the "addr:*" tags. They seem to be a solution for my problem. But in the data I inspected the "addr:state" tag is missing at nearly every point or polygon (at least in the Berlin map I downloaded from Geofabrik).

Is there any (relatively none complex) way or possibility to get these information too from the data without doing geocoding. Or is geocoding the only possibility to solve my problem.

I know that there are mighty tools like Nominatim, that would solve the problem. But I would be happier if there was a more simple way to do that.


texts: Extracting LiDAR metrics with TerraScan?

texts:
I want to extract from a las file statistics such as height percentiles (percentile 95, for example) among others.

I can do it with Fusion (using Cloudmetrics and Gridmetrics tools), but I would like to do it with TerraScan. Is it possible?


texts: Following symbolic links in QGIS to load data on Linux?

texts:
Can QGIS 2.4 follow Linux soft links?

I want to add a number of raster layers from a store of NED30m elevation data, and thought I'd do it by creating a collection of softlinks in one directory to the elevation files nested in the store.

However, when I try to add a layer from one of the links, QGIS 2.4 complains:

  Raster layer: somepath/elev_ned30m_links/n46w105.adf is not a supported
  raster data source

However, n46w105.adf is the result of a linux ln -s command:

  ln -s /Mapstore/ned30m/grid/n46w105_1/w001001.adf

QGIS will open w001001.adf just fine if approached directly, but that means digging down into each tile's directory structure, while the softlinks can be created programmatically with a script.


texts: Add Polygon Attributes to Points error QGIS

texts:
I am trying to run the Add Polygon Attributes to Points algorithm from the Processing Toolbox in QGIS but I run into an error. All my points are either inside the polygon or exactly on the border. I have a feeling that those points exactly on the border might not belong to a polygon. Any suggestions or alternatives?

This is the error message:


texts:
Currently if I export my map to a PDF in ArcMap 10.2, I end up with a PDF that is like 32 MB of data. Our MS Outlook email system caps us off at 15MB per email, so I end up having to open the newly exported PDF map and selecting 'Save As' and 'Reduced size PDF' from the adobe menu.

My question is, is there a way to export a reduced size pdf of my map from ArcMap directly without sacrificing the overall detail and clarity of the image?

My current settings for the export are as follows: 600dpi, Best

I typically export maps with a lot of small print on them and I need the clarity to be fantastic, when I save the 32MB pdf as a reduced size pdf I typically end up with something right around 1MB which is great in terms of size. I just feel like there must be a way to streamline this process.


texts: ArcGIS DissolveManagement Memory Problems

texts:
beside other tools I am using the dissolve tool in a python script. The input feature class (polygon) has about 140 MB. Running the script with python 32 bit (in CMD and in ArcGIS Desktop 10.2.2) fails due to memory problems. Of course the script succeeds with python 64 bit (only in CMD). Well, I made some tests:

If I start the dissolve process manually using the "standard" Dissolve Tool (Generalization Toolbox) in ArcGIS Desktop (using 64 bit Background Geoprocessing), the ArcMap.exe uses 200-300 MB RAM, the RuntimeLocalServer.exe uses about 1 GB RAM. Hence, everything is fine with 32 bit. The process succeeds.
If I integrate the Dissolve Tool in the Modelbuilder and start the model, the RAM that is used by ArcMap.exe is raising up to 1.2 GB, the RuntimeLocalServer.exe seems to be inactive (No background geoprocessing?). However the tool succeeds.
If I start the tool via python script in the command line (or executed in an own Toolbox in ArcGIS Desktop), the ArcMap.exe process is raising up to 2.7 GB RAM and, of course, it fails with python 32 bit.

Why does the same tool with the same input data have a different use of memory in context of ArcMap Toolbox, Modelbuilder and python script via CMD (arcpy)????
Currently I have to run the script with python 64 bit via CMD because the memory is exploding if the script runs within an (own) ArcGIS Desktop Toolbox. Unfortunately I have to change also some other tools (e.g. ET Geowizards) that are not running under 64 bit, too.

Does anybody have an idea why the use of memory is "low" if this tool is only running with 32 bit using the ArcGIS Toolbox?? (or the other way: why is the memory exploding if I am not using the data management/generalization/dissolve tool?


texts: Difference between DynamicMapServiceLayer and TiledMapServiceLayer

texts:
I'm a bit of a noob to Esri still so I have a few questions regarding these two types:

What is the difference between DynamicMapServiceLayers and TiledMapServiceLayers?   Do separate services need to be published to use each type for a single map?  Are both valid types one can use in Esri's FLEX API?


texts: Initalize plugin with current working directory?

texts:
Upon starting Qgis, a (currently under construction) plugin needs a means of checking the current project/legend/registry.

This checksum would ideally occur upon startup automatically.  As the plugin has several buttons which function based on the current project's save location (if one exists).

Example:

Rather then having to click a button in the plugin toolbar every time Qgis opens, the directory resource could be available to the plugin upon initiation and thus establish the credibility of the current project with respect to the plugins functionality.

Important Note:

[The button itself currently has no issue finding the working project directory after startup.  Yet since the qgis-plugins seem to be loaded before the 'last saved' project in Qgis initialization - I have not been able to find a means of automating this process]


texts: driving time to the nearest facility using pgrouting

texts:
I'm using pgrouting 2.0 and postgis in postgresql.
My problem is to calculate the driving time to the closest facility for each node of my OSM derived network.
These are the steps to be implemented (if I'm not mistaken..):
1. Find the closest node of the network for each facility (let't call them "nodes_fac");
2. Calculate iteratively the driving distance to all nodes (or within a given distance) from each of the previous nodes_fac writing the results in a table;
3. Filter the results using the minimum operator to select the driving time to the closest facility (the result should be a table including the nodes geometries with the minimum values);

Step 1: I think that this code can work (posts is the table of facilities):

SELECT posts.id AS post_id,
  (SELECT nodes.id AS node_id
   FROM  ways_vertices_pgr as nodes
   ORDER BY posts.geom <#> nodes.the_geom LIMIT 1)
FROM post_offices as posts
ORDER BY post_id;

Step 2: The following code could be used, but I don't know how to loop it to produce a table with the driving distance for each node_fac:

SELECT ways_vertices_pgr.id,
    ways_vertices_pgr.the_geom,
    basin.cost
   FROM ways_vertices_pgr
     JOIN ( SELECT pgr_drivingdistance.id1,
            pgr_drivingdistance.cost
           FROM pgr_drivingdistance('
      SELECT gid AS id,
          source::integer,
          target::integer,
          avg_cost_min::double precision AS cost
          FROM ways'::text, 22938, 10::double precision, false, false)) basin ON ways_vertices_pgr.id = basin.id1;

avg_cost_min is a calculated column with the driving time for each network segment, 22938 is one of the node_fac, 10 is the maximum driving time.


texts: QGIS 3.16 Could not create print device error (GEOPDF)

texts:
I'm having issues exporting a georeferenced map from QGIS using version 3.16.0-Hannover. I'm getting a "could not create print device" error (more information below).
What I'm trying to do:
I'm trying to export a limited view (approximately 50x50km) of a larger data set (vector data of the whole of Sweden) to a georeferenced PDF in order to view it in map software like "Avenza maps" etc.
My issue is the following:
Using QGIS 3.16 I get an error upon trying to export a geopdf from QGIS layout manager. The error reads "could not create print device".
Link to image of error: https://i.imgur.com/AgIgMEo.png
The export settings I use are the defaults.
Link to image of export settings: https://i.imgur.com/KbPTPyI.png
Example project:
Here's a link to an example project giving me the error (new projects with the same dataset and layout get the same error).
Project link: https://mega.nz/file/S1x0FAYL#SuRhQjFUY_tbyO3iJrLiPsKlEO7G6H4VZ_-_kmwK5nc
Here's what I know/ think so far:

I've tried to google the issue but the results provided there has not been helpful for me.
Exporting a geopdf directly from map view gives me a similary error ("could not save the map to pdf")
Saving a non georeferenced pdf from the same project/ layout and using the same dataset works fine.

What could be the issue?


texts: Iterating through sub-folders within root folder to perform operations on raster files

texts:
I am trying to figure out a way to have my model perform its functions on the rasters contained inside the sub-folders of the root folder. I can't use the default iteration function in model builder because it only looks for rasters in the root folder. This is what my model looks like currently.

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Script arguments
Basins = arcpy.GetParameterAsText(0)
if Basins == '#' or not Basins:
    Basins = "Basins" # provide a default value if unspecified

# Local variables:
Basins_Project2 = "C:\\Users\\straussk\\Documents\\ArcGIS\\Default.gdb\\Basins_Project2"
UnmaskedData = "U:\\SWEModelConstruction\\UnmaskedData"
Raster = "U:\\SWEModelConstruction\\UnmaskedData\\Raster"
SetNull_Rast2 = "C:\\Users\\straussk\\Documents\\ArcGIS\\Default.gdb\\SetNull_Rast2"
Output_Name = "C:\\Users\\straussk\\Documents\\ArcGIS\\Default.gdb\\ZonalSt_Basins_3"
ZonalSt_Basins_3 = Output_Name
ZonalSt_Basins_3__6_ = ZonalSt_Basins_3
Raster_Output = ZonalSt_Basins_3__6_
Basin_Raster_Output = "U:\\SWEModelConstruction\\UnmaskedData\\%Name%.xls"
Name = "Raster"

# Process: Project
arcpy.Project_management(Basins, Basins_Project2, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "", "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]", "NO_PRESERVE_SHAPE", "")

# Process: Iterate Rasters
arcpy.IterateRasters_mb(UnmaskedData, "", "", "NOT_RECURSIVE")

# Process: Set Null
arcpy.gp.SetNull_sa(Raster, Raster, SetNull_Rast2, "\"Value\" = 55537")

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(Basins_Project2, "Drainage_A", SetNull_Rast2, Output_Name, "DATA", "ALL")

# Process: Calculate Field (7)
arcpy.CalculateField_management(Output_Name, "Mean", "[Mean] /1000", "VB", "")

# Process: Calculate Field (8)
arcpy.CalculateField_management(ZonalSt_Basins_3, "Std", "[Std] /1000", "VB", "")

# Process: Calculate Field (9)
arcpy.CalculateField_management(ZonalSt_Basins_3__6_, "Sum", "[Sum] /1000", "VB", "")

# Process: Table To Excel
arcpy.TableToExcel_conversion(Raster_Output, Basin_Raster_Output, "ALIAS", "CODE")

My rasters are placed in folders categorized by their year and date so simply inputting them into one root file is not possible. Is there a way to perform this function on a file system like the following U:\SWEModelConstruction\UnmaskedData\2010\February\15  where the raster files are located on the 15, not the year or month they are just sub-folders to categorize data.


texts: Data needed for wolf habitat/potential distribution in Pacific NW region

texts:
Hello fellow GIS’ers,

I am assisting a nature writer with maps for a book he’s writing on wolves in the Pacific NW (Titled “Wolves in the Land of Salmon”).  I’m writing therefore in the hopes that some of you may know where I can find spatial data on wolves.  My hope is to map wolf subspecies ranges, source populations and regional dispersal corridors, regional historic and current wolf ranges, and projected regional population hubs.  I know of a few habitat suitability studies that have been undertaken (i.e.  Washington Fish&Game, Carrol et al) but so far have been unable to track down much useful GIS data.

I apologize for any cross posting.

Thanks in advance to any leads any of you may have!


texts:
I have downloaded climate forecast data in NetCDF 3.6 format for a series of models and time resolution of one day. The data is one file per model. Eventually, I would like to have an ensemble file with aggregated data over the time axis, e.g. monthly or quarterly. Maybe it gets clearer when I just list some names of the files I have:

rcp85_CNRM_CCLM_RR_daily_mm_2100.nc
rcp85_CNRM_RCA_RR_daily_mm_2100.nc
rcp85_EC-EARTH_CCLM_RR_daily_mm_2100.nc

I expect all the files to have the same extent and resolution; gdalinfo excerpts of arbitrary two files https://pastebin.com/UtxXDDGp support my assumption.
Being a complete noob, I assume that the ensemble would just be the arithmetic average, but I do not know how to average the "payload" axis of rasters over different layers nor over its timeseries/bands.
The manual is not helping me and Ecosia and OS-GEOS search was not helpful either.
Is this possible with QGIS and how would I do it?
About of my QGIS installation gives:
QGIS version                3.14.0-Pi           QGIS code branch           Release 3.14
Compiled against Qt         5.15.0              Running against Qt         5.15.0
Compiled against GDAL/OGR   3.1.2               Running against GDAL/OGR   3.1.2
Compiled against GEOS       3.8.1-CAPI-1.13.3   Running against GEOS       3.8.1-CAPI-1.13.3
Compiled against SQLite     3.32.3              Running against SQLite     3.33.0
PostgreSQL Client Version   12.3                SpatiaLite Version         4.3.0a
QWT Version                 6.1.5               QScintilla2 Version        2.11.3
Compiled against PROJ       7.0.1               Running against PROJ       Rel. 7.0.1, May 1st, 2020
OS Version                  openSUSE Tumbleweed This copy of QGIS writes debugging output.
Active python plugins       processing;
                            MetaSearch;
                            db_manager


texts: I've been informed by my employer of a new policy where I am not to put my name on maps. I've always put my name on my maps in the past and have never run across this policy. Is it uncommon in the industry to add the map author's name to maps?


texts: Importing GeoJson data to postgresql using ogr2ogr command

texts:
I am using following ogr2ogr command to import GeoJson data to Postgres DB. The command runs successfully but the geometry column is null.

ogr2ogr -f PostgreSQL PG:"dbname=test user=postgres password=test" sample3.geojson -nln test.sample3

The geojson file is big so putting some starting lines.

{ "type": "FeatureCollection", "properties": {
"bounds": "0.000001,0.000000,6.343452,51.486215",
"center": "3.779297,51.013722,11",
"description": "sample3.mbtiles",
"format": "pbf",
"json": "{\"vector_layers\": [ { \"id\": \"motorwaygeojson\", \"description\": \"\", \"minzoom\": 11, \"maxzoom\": 11, \"fields\": {} } ],\"tilestats\": {\"layerCount\": 1,\"layers\": [{\"layer\": \"motorwaygeojson\",\"count\": 21144,\"geometry\": \"LineString\",\"attributeCount\": 0,\"attributes\": []}]}}",
"maxzoom": "11",
"minzoom": "11",
"name": "sample3.mbtiles",
"type": "overlay",
"version": "2"
}, "features": [
{ "type": "FeatureCollection", "properties": { "zoom": 11, "x": 1023, "y": 1024 }, "features": [
{ "type": "FeatureCollection", "properties": { "layer": "motorwaygeojson", "version": 2, "extent": 4096 }, "features": [
{ "type": "Feature", "properties": {  }, "geometry": { "type": "LineString", "coordinates": [ [ 0.000000, 0.000043 ], [ 0.000000, 0.000000 ], [ 0.000000, 0.000043 ] ] } }
] }
] }
,
{ "type": "FeatureCollection", "properties": { "zoom": 11, "x": 1024, "y": 1024 }, "features": [
{ "type": "FeatureCollection", "properties": { "layer": "motorwaygeojson", "version": 2, "extent": 4096 }, "features": [
{ "type": "Feature", "properties": {  }, "geometry": { "type": "LineString", "coordinates": [ [ 0.000000, 0.000043 ], [ 0.000000, 0.000000 ], [ 0.000000, 0.000043 ] ] } }
] }
] }

As you can see I have nested FeatureCollections, so what needs to be done to handle this scenario?

My postgres version is 9.3 and gdal version is 2.2.1


texts: What is the simplest way to get the average of multiple raster images while ignore nodata?

texts:
I am looking for the simplest approach to take the average across multiple raster images, which contain a lot of nan values. The mean value should ignore any nan. The gdal_calc does not work, because it outputs nodata when any nodata value is encoutered.

Is there any simple approach with GDAL or python (with any package)? Prefer not to use ArcGIS.


texts:
Using a python script, i would like to read through a table which features 27 columns. For each row i want to store the value of each column in a variable, so that i can pass the values to a subsequent script as input parameter.

here's what i have so far:

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

workspace = r"X:\Auftrag\2015\1415154_Armasuisse_GHK\04_Arbeitsordner\VoF\Armasuisse_Vof.gdb"
arcpy.env.workspace = workspace
arcpy.env.scratchworkspace = workspace
arcpy.env.overwriteOutput = True

Inputdaten_txt = "Batchtable"

Attributes = arcpy.ListFields(Inputdaten_txt)
Attribute_List = []
for field in Attributes:
    Attribute_List.append(field.name)

with arcpy.da.SearchCursor(Inputdaten_txt, Attribute_List) as cursor:
    for field in Attribute_List:
        for row in cursor:
            OBJECTID = ("{0}".format(row[0]))
            Temp = ("{0}".format(row[1]))
            CalculationOption = ("{0}".format(row[2]))
            ElevationModel = ("{0}".format(row[3]))
            DrainageNetwork_rainfall = ("{0}".format(row[4]))
            RelativeDrainageElevation = ("{0}".format(row[5]))
            Elevation_of_input_drainage_network = ("{0}".format(row[6]))
            ConstantDrainatePlus = ("{0}".format(row[7]))
            HydrographFile = ("{0}".format(row[8]))
            CoordinateFile = ("{0}".format(row[9]))
            Roughness = ("{0}".format(row[10]))
            FlowBarrier = ("{0}".format(row[11]))
            DamFailure = ("{0}".format(row[12]))
            ContinueCalculationFrom = ("{0}".format(row[13]))
            ContinueTime = ("{0}".format(row[14]))
            ContinueInundationAsWaterDepth = ("{0}".format(row[15]))
            DurationOfSimulation = ("{0}".format(row[16]))
            SavingIntervall = ("{0}".format(row[17]))
            OutputSwitch = ("{0}".format(row[18]))
            MapUnits = ("{0}".format(row[19]))
            HeightUnits = ("{0}".format(row[20]))
            drainage_units = ("{0}".format(row[21]))
            MaximumExchange = ("{0}".format(row[22]))
            Output = ("{0}".format(row[23]))
            intermed = ("{0}".format(row[24]))
            LicenseCode = ("{0}".format(row[25]))
            Execute = ("{0}".format(row[26]))

            # continue by calling another python script

I'm quite sure that i could be done easier, for instance i would like to automatically generate each variable (e.g. var_1 ...var_27) using a loop that goes through a list of these variables obtained with the searchCursor function. But whenever i tried, i got stuck at the point where my variable ends up being a "SearchCursor object at 0x14A522A0" instead of the field value in the table. Can anyone tell me the easiest way to get this task done?


texts: ArcGIS planarize vs autofixing "must not intersect or touch interior" topology

texts:
trying to fix topology errors in OSM here.  I am wondering whether ArcGIS planarize is essentially the same as creating a topology that says a class must not intersect or touch interior.  Then autofixing all those errors?

One difference I can see is that the topology allows you to use multiple classes; thus you shouldn't break relationships to other classes by a topology fix whereas it might happen in planarize when e.g. two points closer than cluster tolerance are snapped together.

Can anyone enlighten me as to what is going on here?

I think I may have found out already.  ArcMap hangs trying to autofix 50,000 topology errors, it may or may not finish the job but I didn't wait to find out.  It can handle a planarize though.  That suggests the algorithms used are different.


texts: How to get LonLat Points between a LineString or 2 LonLat Points

texts:
i have a marker and i want to move him on a LineString. So i know start LonLat and end LonLat and also the time between.

I tried to calculate it my self with "getLayerPxFromLonLat" and then count x+1 y+1 and after that i switch them back with "getLonLatFromPixel". So i can calculate it my self (e.g. 10 Seconds = 5 Points between start and end) but when i the map moves the function got wrong data.

I hope u have a idea for that Problem or an Idea how to get Points on a LineString.


texts: ArcMap arrow symbol for point: how to make the arrow end starts from the point instead of sitting in the middle of point?

texts: I am trying to use arrows to display points. The arcmap worked nicely to rotate arrow according to the angle field (layer property/symbology/advanced/rotation), however the middle of the arrow sit on the point. I want the arrow to start from point. How can I do it?


texts: I was developed map application using "arcgis runtime for wpf" and this app works fine on my pc, but when I published my project with visual studio and install the application on the client machine the application go down. what is the simplest and the correct way to deploy applications uses Esri controls like 'map and legend'?


texts: Polyline to raster cell size too small to run?

texts:
I'm trying to run a polyline to raster tool and would like the cell size of the outputted raster to match my DEM grid cell size of 1/3(0.0000925925926 degrees) or 1 arc second (0.000277778 degrees).

However the tool is taking forever to run on ArcMap, or fails to run on ArcCatalog.

Is the pixel size just too small for ArcGIS to convert, or is there another way I can achieve this?

I'd like to avoid aggregating the DEM grid.


texts:
Does any body have an idea about why *.tiff images are not displayed when added as raster layer to QGIS after adding a vector layer which has a polygon geometry?

SOLVED !!
In fact the tiff is displayed but not show in the canvas due to a mislocation.
a "(QgisInterface)->zoomToActiveLayer()" after adding the raster layer should bring the image to display ^^
Thanks everybody !


texts:
I created this script to merge all of the polygons in a shapefile into a single polygon. I am using Python 2.7.5, Shapely 1.3.0, and Fiona 1.1.5.

def merge(path_in, path_out):
    #Import Libraries
    import fiona
    from shapely.geometry import mapping, shape, MultiPolygon, Polygon
    from shapely.ops import cascaded_union
    from fiona import collection

    #Import shapefile
    with collection(path_in, "r") as input:
        input_driver = input.driver
        input_crs = input.crs
        input_schema = { 'geometry': 'Polygon', 'properties': { 'INT_DATE': 'int' } }

    #Open the file to be written
        with collection(path_out, "w", driver=input_driver, crs=input_crs, schema=input_schema) as output:
            shapes = []

            #Append all of the polygons to a list
            for f in input:
                shapes.append(shape(f['geometry']))

            #Merge all the polygons
            merged = MultiPolygon(cascaded_union(shapes))

            #Write the new shapefile
            output.write({'geometry': mapping(merged), 'properties':{'INT_DATE':'4'}})

if __name__ == "__main__":
    import sys
    merge(sys.argv[1],sys.argv[2])

When I run this program it returns no error message and outputs a shapefile with the attribute table I expect. However, the polygon is not there. Does anyone have any idea what is causing this? I have relied heavily on the Fiona users manual and "GIS with Python, Shapely, and Fiona" by Tom MacWright.

Is there a better way to use Shapely than using Fiona?


texts: How to get coordinates of label position or make static labels for each feature

texts:
How can I get the coordinates of generated labels on Canvas? Or is there a way to make static labels for each feature? I want the labels to not move when the canvas changes scale or moves.

My situation: my plugin loads the pipelines (LineString) into memory layer. Label for the pipelines must be parallel and above the line. After first rendering of the canvas I need to remember the position of these labels. The ideal variant would be to get the x, y coordinates and angle of rotation of each labels and put them into additional fields of layer.

I'm interested in only stage of extraction of coordinates.


texts: Writing In-memory IFeatureLayer to Shapefile using ArcObjects (Engine SDK)

texts:
I've been working on this for a little while now, and can not seem to get this thing working.

I've managed to export a geodatabase feature class and shapefile to a shapefile, but what I'm struggling with right now is writing an in-memory feature layer (one that is created based on a selection of features) to a shapefile/geodatabase.

If anyone can help me with a code snippet (doesn't matter the language), or some links to some dark corner of the ESRI docs, I would be much obliged.

UPDATE

I used the following code to create the IFeatureLayer:

IFeatureSelection featureSelection = (IFeatureSelection)layer;
if (featureSelection.SelectionSet.Count > 0)
{
    IFeatureLayerDefinition2 featDef = (IFeatureLayerDefinition2)featureSelection;
    IFeatureLayer fLayer = featDef.CreateSelectionLayer(
        String.Format( CultureInfo.InvariantCulture,
            "{0}_Selection", layer.Name ),
        true, "", "" );

    featureSelection.Clear();

    if (fLayer != null)
    {
        _hook.FocusMap.AddLayer( fLayer );
        _hook.ActiveView.Refresh();
    }
}

This creates and adds a new layer to the TOC.

What I'm missing now is taking that Feature selection and saving this to a Shapefile (for now).

The code example I've been working with is from esri code snippets but I'm getting stuck on the following:

IDataset sourceDataset = (IDataset)_layer;

IWorkspace sourceWorkspace;
IPropertySet propertySet = sourceDataset.Workspace.ConnectionProperties;
sourceWorkspace = sourceDataset.Workspace.WorkspaceFactory.Open( propertySet, 0 );

if (sourceDataset.FullName != null)
{
    // HERE - Invalid Cast Exception.
    IWorkspaceName sourceWorkspaceName = (IWorkspaceName)sourceDataset.FullName;
}

The full exception is here:

  Unable to cast COM object of type 'System.__ComObject' to interface
  type 'ESRI.ArcGIS.Geodatabase.IWorkspaceName'. This operation failed
  because the QueryInterface call on the COM component for the interface
  with IID '{FADD975C-E36F-11D1-AA81-00C04FA33A15}' failed due to the
  following error: No such interface supported (Exception from HRESULT:
  0x80004002 (E_NOINTERFACE)).

I guess the exception is thrown because the FullName property though is not empty is not compatible because of the type of COM object it is wrapping.

Any thoughts?


texts: How to extract most probable flow line from multiple lines?

texts:
I have to apply uncertainty principle to computing flow accumulation using Monte Carlo method. I generate raster with random numbers which represents mistakes in creating DEM. Then I apply those rasters on DEM in iterations. For each DEM with added random numbers raster I compute flow accumulation. In each iteration it generates slightly different lines of flow accumulation. At the end of process I merge these lines into one raster. It looks like this: multiple flow lines

What I have to do is extract most probable flow line from this data. I use ArcGIS10 but it does not matter which software should be used. In ArcGIS I tried to use thin function or extract by attributes which gives me pretty good results. I also converted every line into polyline, merged them and applied buffer. Then I used thin function on buffer. This gives me also very good result.

I want to ask if you can help me with finding another solutions for extracting most probable line?


texts:
Where can I put the 'crossdomain.xml' file in the MapProxy virtual enviroment ?

[UPDATE] [NEED TAGS for question: mod-wsgi, ogcserver too]
Need to be more specific:

i'm using mapproxy 1.4 in mod-wsgi.
The same problem I had found in OGCServer (mod-wsgi).
Where have i to put the crossdomain.xml in a mod-wsgi server like Mapproxy or OGCServer?

Thanks again.


texts:
I'm entirely new to QGIS and mapping services in general.

I add my data via CSV to create a new map layer (following the QGIS heatmap tutorial), I select Longitude in my X coordinate, and Latitude in Y, as shown below. If I don't do that, my data does not appear as it does when geocoded onto map services (such as BatchGeo). This seems counter-intuitive to me.

From there, I am attempting to add a basemap underneath my data points (ideally underneath a heatmap I've made - not sure how to do this either) but I seem to be having trouble. I have attempted using the OpenLayers PlugIn instructions found here (http://maps.cga.harvard.edu/qgis/wkshop/basemap.php), but I keep getting an error when trying to add a Google or Bing Map. The message says "OpenLayers Plugin: Printing and rotating of Javascript API based layers is currently not supported!"

I have since found that QuickMap Services seems to be a preferred plugin, but I seem to have limited map options, and when i do pull in a map, my data points don't seem to show up where they need to be (looks like lat/long is reversed, which might relate back to how I import the data?).

Does anyone have any solutions for beginners for how to add a basemap or even alternate tutorials for creating heatmaps (on top of a basemap)?


texts: Superimpose grid on points, then aggregate points within each grid cell

texts:
I have a vector of presence/absence values (1,0) and a matrix of x,y coordinate data (projected) specifying the spatial position of each value (about 25,000 in total). I'd like to aggregate the point data into 'cells' and calculate percentage values (% of presence) for each cell. However, I want the cells to be arbitrarily positioned with respect to the point data (i.e., not determined by census tract or any other factor), and with a fixed uniform size (0.5 square kilometers each).

So essentially, I'd like to superimpose a fixed-size grid over my point distribution and then calculate the percentage of 1s (presence) in each cell of the grid. I then want to repeat the proceedure many times, each time randomly altering the position of the grid with respect to the points, and recalculating the percentages of each cell. Ultimately, I will perform localized spatial tests on the percentage data.

I realize there are probably multiple ways to approach this using many tools. I have some familiarity with R and the 'spdep' package and basic famililarity with ArcGIS. I'd be grateful if someone could at least point me in the right direction for R or ArcGIS functions that may be useful.


texts: Query against FeatureClass for layer using ArcObjects?

texts:
This is some code that I am using to try to get a working IFeatureLayer from a subset of data from a database table. The fourth line containing "fcWkspc.OpenFeatureQuery" throws and error complaining about an "The operation was attempted on an empty geometry.".

Should also mention that this code was lifted directly from ESRI documentation example code at the following link:

http://resources.esri.com/help/9.3/arcgisdesktop/arcobjects/esrigeodatabase/IQueryDef_Example.htm

Can anyone tell me what I might be doing wrong?

// m_wkSpace2 is of type IWorkspace2
IFeatureWorkspace fcWkspc = (IFeatureWorkspace)m_wkSpace2;

IQueryDef qDef = fcWkspc.CreateQueryDef();

Def.Tables = "dbo.TABLE_NAME";
qDef.WhereClause = "Name1 = '" + var1 + "' AND Name2 = '" + var2 + "'";

IFeatureDataset featureDataset = fcWkspc.OpenFeatureQuery("MyJoin", qDef);
IFeatureClassContainer featureClassContainer = (IFeatureClassContainer)featureDataset;
IFeatureClass featureClass = featureClassContainer.get_ClassByName("MyJoin");

IFeatureLayer layerForMap = (IFeatureLayer)featureClass;


texts: Using CSS var function in GeoServer CSS styling

texts:
I wanted to check if I could use something like the CSS var() function with the GeoServer CSS plugin styles.
Something to look like:

    :root {
        --land-color: #f2efe9;
      }
    /* @title Land */
    [featurecla = 'Land'] {
        background-color: var(land-color);
      }

Except obviously for more than one variable.
I came across this question but could not get the only suggested solution there to work.
This is what I tried:

/* &env=land-color:#dd0000 */
/* @title Land */
[featurecla = 'Land'] {
fill: env('land-color');
stroke: #f2efe9;
}

EDIT

The result I want to get is setting a color (or any other style variable, like halo-radius or font-size) once and the use it by name in entire stylesheet, making it easier to edit it across the document.

I keep getting a default white polygon with black stroke when using this. The CSS validates properly but does nothing.


texts: How to weigh the size of a spatial event with respect to the number of its hotspots?

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

My problem: I want to find in the map the regions (or islands) which values exceed a given threshold, e.g. 0.3 (we can call them the hotspots). I then want to calculate the number of pixels or cells that fall within the boundaries of that island. Given that there might be more than one island, I could end up with:

Situation 1: 1 island with, say, 1000 pixels or cells;
Situation 2: 2 islands with 500 pixels or cells each;
Additional combinations.

Since I must be able to distinguish between the various situations, I was thinking of somehow weighing the total number of cells with respect to the number of islands. How could I do this, other than considering the average number of cells (which I see as a trivial solution)?

I am asking this since if I were to evaluate the impact of those hotspots on a system that is embedded in the map, there might be huge differences between having one big hotspot located in the center of the map and having multiple, although smaller, hotspots scattered all over the place.

Below is a visual example of what I mean by islands or regions.


texts:
Taking some code from this post, I've been trying to create an 'editable' view in PostGIS.

I'd like to store label placement properties in a separate table, build a view of the data I want to use for my map, and use QGIS to move the labels. The idea is when the labels are updated, the placement coordinate would go into the LABEL_X and LABEL_Y fields - the geometry of the underlying points and several other attributes (in the view, coming from the source tables) would remain unchanged.

However, when all is in place, QGIS throws the following error when trying to save the edits to the view (the edit session can start, labels moved, but cannot save):

Could not commit changes to layer schools_district_map

Errors: ERROR: 1 attribute value change(s) not applied.

Provider errors: PostGIS error while changing attributes:

ERROR:  infinite recursion detected in rules for relation "schools_district_map"

Here is the definition of the view:

CREATE OR REPLACE VIEW public.schools_district_map AS
 SELECT sch.schnum,
    sch.oid,
    sch.abbreviation,
    sch.school_level,
    sch.geom,
    l.label_x,
    l.label_y
   FROM temp_schools_label sch
     LEFT JOIN district_map_labels l ON sch.schnum = l.schnum;

And here are the rules I've applied to make the view 'editable':

--delete rule
create or replace rule "delete_label" as
on delete to schools_district_map do instead
delete from schools_district_map
where oid = old.oid;

--insert rule

create or replace rule "insert_label" as
on insert to schools_district_map do instead
insert into schools_district_map (label_x, label_y)
values (new.label_x, new.label_y);

--update rule
create or replace rule "labels_update" as
on UPDATE TO schools_district_map do instead
update schools_district_map set
label_x = new.label_x
, label_y = new.label_y
where oid = new.oid;

QGIS is then set to display the labels using the label_x and label_y field.

If anyone can spot where I might have left something out, or if there is a glaring oversight on my part, let me know.


texts:
Essentially I want to union multiple (20) congruent line shapefiles, where overlapping line segments in the input shapefiles are combined into a single resulting output shapefile. But I can't union lines. So I'm looking for an alternative. The desired result, like union, would break all the overlapping lines into separate features and preserve the input attributes.

Merge doesn't work for me because the inputs are preserved as separate (overlapping) features in the resulting shapefile. I tried dissolving this result, as well as Split Line at Vertices without getting the desired result.


texts: Need a location based social network data set

texts:
I'm looking for a data-set of Geosocial network consisting of GPS trajectories, check-in, geotag contents and stuff like these available for research purposes.

I'd really appreciate it if you can point me towards such a reliable data-set.


texts: pgRouting - calculating route from point to point

texts:
I have three shapefiles:

roads(lines),
logistic center(one point feature) and
places(a lot of point features).

I have tried to calculate route from logistic center to places by ArcMap network analyst, but it was pain and I do not successfully made it.

So I installed PostGIS/Postgresql with pgRouting extension. Now I have those shapefiles imported to dabatabase with SPIT tool from QGIS. In database they are shown as 3 tables.
Also I successfully created topology for roads by tutorial found in pgRouting workshop documentation.

For now, I can't find tutorial for my situation, as I want to calculate shortest routes from logistic center to places using roads layer. I can do it by hand, one by one, on QGIS or ArcMap, using their tools, but I have 1500 places, so it would be pain to do it.

My questions:

Is it possible to calculate routes from one point to many with pgRouting (I need only distances for shortest route) using logistic center as source and places as targets?

Is it possible to create a table with information containing target ID and route distance and create a shapefile with geometry?

How to do that?


texts: style manipulation using Google Earth api

texts: I am creating data in QGIS and then serving it using GeoServer and displaying using the Google Earth api. Is there any way to access the styles using a JavaScript api and alter the look and feel of my data, (i.e. changing colors, etc.)?


texts: How to export multiple MXD's using python?

texts:
I have several mxd's that use data driven pages to create an updated mapbook every month. 15-18 mapbooks and 163 pages each. This usually takes me two days and I'd like to automate this. I can export one with python, but more than one is stumping me.
I export one using this in Python:

map1 = arcpy.mapping.MapDocument("map1 mxd location")

map1.dataDrivenPages.exportToPDF(r"folder where I want pdf stored", "RANGE", "1-163","PDF_MULTIPLE_FILES_PAGE_NAME")

del map1

It works great for one map, but how can I export map1, map2, map3 etc.. to separate folders with one executed command?

Just to be clear, I know I need to specify each mxd and export location individually. What I'm looking to learn is how to enter all that information at once, double tap my enter button and walk away.


texts: Digitizing Airport approach-zone with z information in QGIS?

texts:
I'm trying to create a vector layer in QGIS out of this Z-dimensioned approach-zone for an ongoing multi-criteria analysis:

* The oval shaped base should be 59,56 meters high (Z)

For every "ring" z-value should be increased by 10 meters, ending
at 159,56 meters (Z)
The cone shaped things are flat plane 163/164 meters up until the 130
meter mark, then drops until it reaches 59.56.

Could anyone with more experience with QGIS-digitizing with z-values point out a direction?

What painting tools to use? Right now I'm manually painting them, but results looks horrible and edgy.

note: I've already geo-referenced the above image.


texts: Spatial autocorrelation coefficients interpretations

texts:
I have this result for SDM estimates in R with R2 of 54%. Is the lambda (= -2.04) which is the spatial autocorrelation coefficient a plausible value? I have read in a book that [Elhorst, 2014, P.14], the spatial autocorrelation coefficient should be bounded between -1 and 1.

    > summary(ERVM)
Spatial panel fixed effects lag model

Call:
spml(formula = fmerv, data = spa.sakk, index = c("id", "months"),
    listw = ERV, model = "within", effect = "time", lag = TRUE,
    spatial.error = "none")

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max.
-0.785173 -0.018300 -0.000716  0.016336  1.112927

Spatial autoregressive coefficient:
       Estimate Std. Error t-value            Pr(>|t|)
lambda   -2.403      0.111   -21.6 <0.0000000000000002 ***

Coefficients:
                Estimate Std. Error t-value             Pr(>|t|)
spa.sakk$credk   0.00298    0.04250    0.07              0.94404
spa.sakk$chexk   0.04923    0.03896    1.26              0.20633
spa.sakk$uninfk  0.23803    0.26140    0.91              0.36250
spa.sakk$gdpk    0.17135    0.33038    0.52              0.60401
spa.sakk$chintk -4.53205    1.10230   -4.11             0.000039 ***
spa.sakk$lagret -0.26580    0.01386  -19.18 < 0.0000000000000002 ***
slervcredT       0.04790    0.14141    0.34              0.73480
slervchecxT     -0.32075    0.12271   -2.61              0.00895 **
slervuninfT      7.46349    4.35165    1.72              0.08633 .
slervgdpT        1.11652    1.01416    1.10              0.27093
slervchintT     -6.41103    1.81507   -3.53              0.00041 ***
sllagret        -0.04978    0.05205   -0.96              0.33882
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


texts: How to modify a DTM in order to show the sea level around an island?

texts: In ArcGIS 10.1 (advanced license), I have a DTM representing a small island. Of course, off-coast cells are given NoData. Now, I would like to have those off-coast areas set to 0, i.e. sea level. I used raster calculator (with a conditional statement) to set the NoData cells to 0, leaving unchanged the land cells. The problem with that is that the sea level is limited to the extent of the raster, but I would like to have a larger portion of sea "represented" (i.e., well outside the raster's extent).


texts: Creating a symbology value for "all other values" in QGIS 3.0 - clarifying question

texts:
This question centers the Symbology tab in QGIS 3.0.

I'm working with categorical data describing event venues in the Chicago area. This is a point dataset which will be overlayed on a street basemap.

Bars
Nighclubs
Theatres
Parks
Restaurants
Galleries
Other

I would like to use the Symbology tool in the Properties window to display "Bars" specifically in red, and "all other values" (nightclubs, theaters, other, etc.) in yellow.

The attribute field I am using is called "Type"

Fundamentally, I have solved this problem using the following process:

Symbology
Column Drop Menu -> Select "Categorized"
"Formula" button
"Type"  =  'Bars' (syntax)
Select "Classify" button
Outcome: classifcation of "0" represented in Yellow, and "1" in Red.
"1" = "Bars", "0" = "all other values/anything that is not a bar"
Change Legend entry for "1" value to read "Bars"
Change Legend entry for "0" value to read "All Other Venues"

Indeed, this process solves the problem of isolating the desired value and representing all other values.

I would however like to clarify... is this the only way to isolate one value and represent all others?  Many of the folks I work with are coming from the ArcMap 9 and 10 environments, which make this process super easy using selectable options in the GUI rather than formulas.

Have I missed a GUI method of completing this task while solving this problem or should I start working with formulas?


texts: ArcPy: Error 000732 feature layer / table view does not exist or is not supported

texts:
I am trying to update a feature class by joining it to a table after deleting columns in the target layer with the same names as columns in the join table. I thought the only way to do this in arcpy was to make the featureclass a feature layer and make the table a table view. I tried my script with out using featur layer or table view and got an error that the datatype was not a raster or mosiac. I then made the feaute layer and table view but now am getting:

ERROR 000732 Dataset does not exist or is not supported.

Here is my script so far:
# script to update Signal dashboard

import arcpy
from arcpy import env

# set up env DB to overwrite outputs
env.workspace = r"C:\Users\Master\Documents\ArcGIS\Projects\Signal_Mgmt\Signal_Related_Apps.gdb"
env.overwriteOutput = True

# set up paths vars, etc
outPath = r"C:\Users\Master\Documents\ArcGIS\Projects\Signal_Mgmt\Signal_Related_Apps.gdb"

# the table used for updating featureclass
table = r"C:\Users\Master\OneDrive\Desktop\Temp out\EOC\Signal_MGMT_Census\AC_zips_rates.csv"

joinTable = arcpy.TableToTable_conversion(table, outPath, "AC_zips_rates")

# Current featureclass published to AGOL
inFC = r"C:\Users\Master\Documents\ArcGIS\Projects\Signal_Mgmt\Signal_Related_Apps.gdb\Zips"

# Get rid of fields from previous update
delFields = ["Count", "Jul_2020_Pop", "Rate_per_1000", "Popup", "Legend", "Bin"]

arcpy.management.DeleteField(inFC, delFields)

# make layer and table views to perform join ops
inField = arcpy.MakeFeatureLayer_management(inFC, "target_lyr")

joinView = arcpy.MakeTableView_management(joinTable, "join_view")

inTab = arcpy.JoinField_management(inField, "Zip_Number", joinView, "Zip")

# export joined tables as new FC
joinOut = arcpy.FeatureClassToFeatureClass_conversion(inTab, outPath, "Signal_Update")

# Remove join from current featureclass so join does not interfre with overwrite
arcpy.RemoveJoin_management(inField)

# Overwrite current featureclass with new featueclass
upFC = arcpy.FeatureClassToFeatureClass_conversion(joinOut, outPath, "Zips")

# delete join field from new featureclass
arcpy.management.DeleteField(upFC, "Zip")

# apply symbology layer to new featireclass
arcpy.ApplySymbologyFromLayer_management(upFC, r"C:\Users\Master\OneDrive\Desktop\Temp out\EOC\Signal_MGMT_Census\Signal_Zips.lyrx")

Am I going about this the wrong way? I can't seem to do a join op either way. I'm using ArcGIS Pro 2.5.


texts: Deleting points based on number of neighbours using ArcGIS?

texts:
I have one point dataset of around 270,000 points. I have another larger point dataset of about 3.5 million points.

I would like to make a 5 km radius around each 270,000 points of the first dataset. If this radius contains less than 3 points (both point datasets together), I would like to delete all points within this radius using ArcGIS.

How can I do this?


texts: I need to show the area of both the water and cropland as comparison, but it gives me only the chart of water

texts:
i calculated the area for two classes with the below code and want to see the bar chart of area for both the classes. When I write the code it gives only the bar chart of one class not as a comparison. need your help
//****************************Area Calculation******************************
//Water Area

var Water = classified.select('classification').eq(0);

//Calculate the pixel area in square kilometer

var area_Water = Water.multiply(ee.Image.pixelArea()).divide(1000*1000);

//Reducing the statistics for your study area

var stat1 = area_Water.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e9
});

//Get the sq km area for Water

print ('Water Area (in sq.km)', stat1);

//Cropland Area

var CropLand = classified.select('classification').eq(1);

//Calculate the pixel area in square kilometer

var area_CropLand = CropLand.multiply(ee.Image.pixelArea()).divide(1000*1000);

//Reducing the statistics for your study area

var stat2 = area_CropLand.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e9
});

//Get the sq km area for Cropland

print ('Cropland Area (in sq.km)', stat2);

var chart = ui.Chart.feature.byProperty(
  ee.Feature(null, stat1)
)

print(chart)


texts: QGIS on docker container: Could not connect to any X display

texts:
I am trying and failing to establish functional environment on docker container to run some QGIS related tasks.

This is the idea: using flask app user enters input data. That information travels to server and it is user to query spatial database. Next, that data is exported to geopackage format and stored to server. After that I want to create qgs project and add geopackage to project, add qml or sld as a style, save it, zip whole folder (with QGIS project and geopackage) and send it to user.

I used docker to do this. I installed qgis/qgis image from docker hub.

The problem starts where I try to instantiate QgsApplication and it says next:

QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-root'
qt.qpa.screen: QXcbConnection: Could not connect to display
Could not connect to any X display.**

This is the way I tried to solve problem but error repeats (also saw on gis.stackexchange):

@app.route('/')
def index():
    vers = Qgis.QGIS_VERSION
    # I got value for vers variable
    ogr2ogr.main(["", "-f", "GPKG", "export.gpkg", "PG:dbname='dbname' host='host' port=5432
    user='username' " "password='password'  tables=schema.table", "-skipfailures", "-
    overwrite", "-progress"])
    gui_flag = False
    app = QgsApplication([b''], gui_flag)
    app.setPrefixPath('/usr')
    app.initQgis()
    if len(QgsProviderRegistry.instance().providerList()) == 0:
        raise RuntimeError('No data providers available.')

    strProjectName = "my_project.qgs"
    project = QgsProject.instance()
    vlayer = QgsVectorLayer("export.gpkg", "layername", "ogr")

    if not vlayer.isValid():
        print("Layer failed to load!")
    else:
        project.addMapLayer(vlayer)

    project.write(strProjectName)
    return render_template('ind.html', vers =vers)

This is Dockerfile I used to build image.
docker image build -t qgis_proj

app.py is a simple flask app with one endpoint as mentioned above.

FROM ubuntu:bionic
FROM qgis/qgis
RUN apt-get update -y
RUN apt-get install python3-pip
RUN apt-get install xvfb xserver-xephyr vnc4server -y
RUN pip3 install pyvirtualdisplay
RUN apt-get update -y
RUN Xvfb :99 -ac -noreset &
RUN export DISPLAY=:99
RUN mkdir /app
WORKDIR /app
COPY requirements.txt /app
RUN pip3 install -r requirements.txt
RUN rm requirements.txt
RUN apt-get install wget gpg nano -y

EXPOSE 5002

ENTRYPOINT ["gunicorn", "-b", "0.0.0.0:5002", "-w", "9", "app:app"]

After building an image an run container
docker container run -d -p 5002:5002 -v /home/cug/download:/app --name download_service qgis_proj


texts: I have to draw all the car park spaces for several streets, I am asking for some advise, how is the best way to draw all those car park spaces? and (of 5 meters) applying some restrictions (Ej. like double yellow line). I am using QGIS and I am not going to have pre-made lines or anything, all just from base maps. (the bay should be 5 meteres each with some areas with the restrictions mentioned before).


texts: Getting ArcMap reference without launching application?

texts:
I need to get hold of ArcMap properties in order to change the data source of the layers within it.

I have developed a similar utility earlier utilizing ICommand, but I was curious if I could get a reference to ArcMap application outside of it and change the layer sources and display the new map.

I am using the following code but it instantly launches ArcMap and after that I have no hold on the application i.e. if I terminate my application, ArcMap is still running in background.

Is it possible to get a silent, invisible reference to ArcMap?

IApplication pApp;
IDocument pDoc;
IMxDocument pMxDoc;
string MXDPath = "C:\\NORMALTEMPLATE\\OLD_Normal.mxd";

pDoc = new MxDocumentClass(); //ArcMap launches here
pApp = pDoc.Parent;

pApp.OpenDocument(MXDPath);
MessageBox.Show("Doc opened");


texts:
I am taking soil samples on an agricultural plot and will be using the results as input to create a semivariogram for kriging in ArcGIS. This plot features raised planting beds in which the test crop is planted. Ideally I am looking to take soil samples using a stratified random approach; with an equal number of samples collected from different management zones within the plot.

I am wondering whether taking samples from the planting beds only will result in valid output from the kriging process. Clearly, using this method, the sample points will be biased in their distribution; linearly, along the planting beds.

Will taking samples from the planting beds only result in a valid interpolation surface?


texts: Load Shapefiles into a remote PostGIS db using command line?

texts:
Can we import a Shapefile from a system which doesn't have PostgreSQL installed (client machine) into a system having PostGIS installed?

I am getting problem with 'shp2pgsql' and 'psql' since both are not recognized on the client machine.

I am trying this using console.


texts: Using R Shiny to create a plot with a shapefile and a raster

texts:
I just wanted to upload a shapefile and a raster in a RShiny app and then execute a script. What I have wrote is the following;

library(shiny)
library(raster)
library(rgdal)
ui=shinyUI(fluidPage(pageWithSidebar(
  headerPanel("Header1"),
  sidebarPanel(
   fileInput('layer', 'Choose Layer', multiple=FALSE, accept='tif'),
   fileInput('shapefile', 'Input Shapefile/AOI',accept=c('.shp','.dbf','.sbn','.sbx','.shx',".prj"), multiple=TRUE)

  mainPanel(
    plotOutput("mapPlot")
  )
)))

server = shinyServer(function(input,output){

  inFile <- reactive({input$layer})
  data <- brick(inFile())
  inFile2 <- reactive({input$shp})
  data2 <- readOGR(inFile2())

  output$mapPlot<-renderPlot({plot(data)
    plot(data2, add=T)})
})
shinyApp(ui = ui, server = server)

What I am getting is the following error

  Warning: Error in .getReactiveEnvironment()$currentContext: Operation not allowed without an active reactive context. (You tried to do something that can only be done from inside a reactive expression or observer.)
  Stack trace (innermost first):

  45: .getReactiveEnvironment()$currentContext
     44: .dependents$register
     43: inFile
     42: brick
     41: server [~/.active-rstudio-document#18]
      5:
      4: do.call
      3: print.shiny.appobj
      2: print
      1: source

I have no idea what I did wrong and I can't find a source that explains how to solve this. Can anyone help me debug it or even help me understand what I did wrong?


texts:
I'm trying to run Whitebox GAT's Import GeoTIFF tool.  The following message pops up each time:

ArcMap's source info says the Tiff I'm trying to import is not compressed or 48/64-bit RGB.  It is a large 32-bit floating point file (~20GB), georeferenced (ie has TFW file) and projected into a StatePlane Coordinate System.  Perhaps the size or projection is causing issues?  I looked at the WhiteboxLog0_0.xml file.  Honestly, it's hard for me to decipher that file in a browser, but the first line is:

"2018-08-30T13:35:07 1535650507206 65 whiteboxgis SEVERE whiteboxgis.WhiteboxGui logException 18 Whitebox 3.4.0. Error Error in MapRenderer"


texts: How to skip drag event to stop a feature being dragged

texts:
I am  new to OpenLayers and I'm trying to drag markers and skip some markers from being dragged, I've tried dragFeature.deactivate(); it will deactivate entire features.

Please check this code :

dragFeature = new OpenLayers.Control.DragFeature(vectorLayer, {
    autoActivate: true,
    onEnter: function (feature) {
        //  example assign : x=1;
        //  if(x==0){ continue draging of a feature }
        //  else{  stop a feature from being dragged }
    },
    onLeave: function (feature) {

    },
    onStart: function (feature, pixel) {
    },
    onComplete: function (feature, pixel) {
    }
});
map.addControl(dragFeature);
dragFeature.activate();

Please help me.


texts:
I have this code:

def main():

    qgisApp = QgsApplication([], True)
    QgsApplication.setPrefixPath('/usr', True)
    QgsApplication.initQgis()

    layers = []

    fileName = "."
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
path = fileInfo.filePath()
rlayer = QgsRasterLayer(path, baseName)

    layers.append(rlayer)

    QgsMapLayerRegistry.instance().addMapLayers(layers)

    img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)

    p = QPainter()
    p.begin(img)
    p.setRenderHint(QPainter.Antialiasing)

    layers = QgsMapLayerRegistry.instance().mapLayers()
    lst = layers.keys()
    renderer = QgsMapRenderer()
    renderer.setLayerSet(lst)

    extent = renderer.fullExtent()
    rect = QgsRectangle(extent)
    # print rect.toString()
    rect.scale(1.1)
    renderer.setExtent(rect)

    # set output size
    renderer.setOutputSize(img.size(), img.logicalDpiX())

    # do the rendering
    renderer.render(p)
    p.end()

    img.save('test.png', 'png')

    QgsApplication.exitQgis()
main()

This only works for pictures that are not 32 bit TIFFs. Is this due to me having an old version of QGIS? How can I get past this?


texts: Nearest Neighbor Based on Network of Lines in QGIS

texts:
Is it possible to incorporate a network into the MMQGIS hub distance analysis? Essentially I have a point shapefile with zip code centroids and I'd like to determine the closest shipping hub (also point shapefile) for each assuming the path must follow a road (a line shapefile).

Unlike creating a road distance matrix I would like to have a single hub identified for each zip code centroid. A further distinction is that neither the zip code centroids nor the hubs actually lie directly on the roads so I'd like to have the source (zip centroids) and destination (hubs) be the closest point on the road/line.


texts: Deleting polygon within polygon of same feature from raster conversion using ArcGIS Pro

texts:
I am currently analyzing some coral reef data. I have PNG files that I have inserted into ArcGIS Pro that look like this zoomed in. The red represents actual coral material while the black is substrate.

I have used raster to polygon to create a polygon feature of this coral that turns out like this:

Instead of creating a polygon layer that has orange polygons that represent only coral material, the raster to polygon tool instead fills in the space inside of coral material polygons that is supposed to be substrate to also be representative of coral material. This is incorrect.
Is there anyway that I can either:

convert the raster to a polygon without filling in these spaces; or
be able to delete these polygons within the actual coral material polygons using a tool or selection feature all though they are the same layer.

I have previously just used editor mode to go through and delete them by hand, by this takes several days as there are so many.


texts: Different results in metrics calculated with LAStools and lidR

texts:
I have canopy height model calculated from high density TLS data on 60 by 200 meter plot.
I tried to calculate voxels with LAStools and lidR and got significantly different results. I was wondering if someone can make clear what is happening.
Lastools script that I used:
lasvoxel -i infile.laz -drop_class 2 -step 0.5 -o outfile.las

number of voxels: 189077

lidR code:
las = readLAS("infile.laz", select = "xyzc", filter = "-drop_class 2 -drop_z_below 0 ")
voxels <- voxelize_points(las, res = 0.5)

number of voxels: 196257

Also calculated canopy height skewness and kurtosis:
LAStools:
lascanopy -i infiles\*.laz -kur -ske -height_cutoff 1.3 -files_are_plots -names -o outfile.csv

Result:
plots        ske       kur
72a-4.laz   1.0905  5.58125
11a-4.laz   0.362   2.594
34-2.laz    0.1675  2.00875
63a-1.laz   -0.3115 2.36

lidR:
library(e1071)

files <- list.files(path= "/files", pattern= "*.laz", full.names = TRUE, recursive = FALSE)
O = lapply(files, function(x) {

  las <- readLAS(x, select = "xyzc")

  z <- las$Z
  z_canopy <- z[z>=1.3]
  skew <- skewness(z_canopy)
  kur <- kurtosis(z_canopy)

  return(data.frame(file=x, skewH = skew, kurH=kur))
})

Result:
plots          ske        kur
72a-4.laz   1.090595768  2.58132381
11a-4.laz   0.362007296  -0.40599745
34-2.laz    0.167542141  -0.991227478
63a-1.laz   -0.311523396 -0.640029907

As we can see, the results for skewness are the same, but kurtosis values are very different.
Can someone please help me understand why there is such a big difference?


texts: Add makeFeatureLayer layer to map in ArcGIS Pro using ArcPy

texts:
I am making a feature layer and applying an expression to essentially apply a definition query to the source feature class.  I then want to add the created feature layer to the map. I am aware of addDataFromPath but I don't want to add the entire feature class. I want to add the layer with a definition query applied. I could save the new layer to file and add it that way but it seems redundant to have the feature class as well as its subset stored on file. Does anyone know how I can add this type of layer to the map in ArcGIS Pro?

c_project=arcpy.mp.ArcGISProject("current")
c_map=c_project.listMaps()[0]
arcpy.MakeFeatureLayer_management ("pts","pts_layer", 'Solution = 1')


texts: Any pgrouting enabled hosting providers?

texts:
I am really interested if any of you know about a hosting provider that offers postgis & pgrouting pre-installed?

If not, can I just go ahead and purchase a linode hosting account and do all the installations myself?

The situation is that I am seeking for a hosting account to launch my pgrouting based application.

Please consider that I am a WIN 7 user with no linux and ubuntu experience.


texts: Access PostGIS Raster layer through mapserver

texts:
I am trying to use a cgi script calling mapserver to display a raster layer that is stored in a PostGIS database. However, only a blank white image is returned and I don't know where the problem is.

The request seems to go through, but here is an error in the postgres log that I don't understand and the result is a white image.

My questions are:

What does the error from the postgres server log mean (see below)?
How can I make sure the white image is not due to a misspecification of the map extent?
Is there something obvious that I am doing wrong?

Here is some more detailed info for this problem.

The mapserver request is:

http://someurl.com/mapserv.cgi?map=../../usr/maps/test.map&mode=map&layers=kahoo&

The error from the postgres server log:

ERROR:  column "o_column" does not exist at character 39
STATEMENT:  select o_table_name, overview_factor, o_column, o_table_schema from raster_overviews where r_table_schema = 'public' and r_table_name = 'kahoo' and r_column = 'rast'

From the mapserver log:

.632473 msDrawMap(): rendering using outputformat named png (AGG/PNG).
.632541 msDrawMap(): WMS/WFS set-up and query, 0.000s
.633883 msDrawRasterLayerLow(kahoo): entering.
.859561 msDrawMap(): Layer 0 (kahoo), 0.227s
.859588 msDrawMap(): Drawing Label Cache, 0.000s
.859594 msDrawMap() total time: 0.227s
.862701 msSaveImage(stdout) total time: 0.003s
.862751 mapserv request processing time (loadmap not incl.): 0.231s
.862759 msFreeMap(): freeing map at 0x266bc60.

The mapfile used is:

MAP
 CONFIG "MS_ERRORFILE" "../../tmp/ms_error.txt"
 DEBUG 5

 NAME "test"
 SIZE 218 152
 EXTENT 739335 2268075 7761335 2284075
 UNITS METERS
 STATUS ON
 PROJECTION
  "init=epsg:26704"
 END

 WEB
  METADATA
    "wms_title"           "WMS Demo Server"
    "wms_onlineresource"  "http://someurl.com/mapserv.cgi?map=../../usr/maps/test.map&"
    "wms_enable_reqest"   "*"
  END
 END

 LAYER
    NAME kahoo
    TYPE RASTER
    STATUS ON
    DATA "PG:host=*** port=*** dbname='***' user='***' password='***' schema='public' table='kahoo' mode='2'"
    PROCESSING "NODATA=0"
    PROJECTION
       "init=epsg:26704"
    END
 END
END

Here are the versions I am using:
mapserver-6.0.3, postgresql-9.1.6, postgis-2.0.1


texts:
I am trying to show a line on map using OpenLayers 3. Here is my code so far:

<script type="text/javascript">
    var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Tile({
                source: new ol.source.XYZ({
                    url: 'http://{a-d}.freemap.sk/T/{z}/{x}/{y}.png'
                })
            })
        ],
        view: new ol.View({
            center: ol.proj.fromLonLat([20.4511,48.91098]),
            zoom: 12
        })
    });
    var points = new Array(
        new ol.geom.Point(20.4511, 48.91098),
        new ol.geom.Point(21.4511, 46.91098)
    )
    var feature = new ol.Feature({
        geometry: new ol.geom.LineString(points)
    });

    feature.getGeometry().transform('EPSG:4326', 'EPSG:3857');
    var vectorSource= new ol.source.Vector({
        features: [feature ]
    });
    var vectorLayer = new ol.layer.Vector({
        source: vectorSource
    });
    map.addLayer(vectorLayer);
</script>

However I cannot get it working. I only see my map which is loaded from source but I do not see any vector line. I have set strokeWidth to 50 to be sure that I will see it but it is not working. I think that there may be some problem with projection.

Any ideas how to get this working?


texts: Unable to connect WFS datasource with ogr, but works with QGIS

texts:
I am accessing some WFS layers via QGIS.
I want to download those same layers via ogr. This doesn't work and I can't find why.

The WFS URL is https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wfs
The layer I want to acess from that WFS is "anb_insp:PS.ProtectedSite-birds"

In the windows command prompt, I put:
ogrinfo WFS:"https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wfs&version=2.0.0&request=GetFeature&typeName=anb_insp:PS.ProtectedSite-birds&srsName=EPSG:31370"

I get "unable to open datasource ... with the following drivers.

I cannot understand why I can access this layer via QGIS, but not through ogr? Am I missing somehting?

PS 1
I am able to open WFS layers from other sources
PS 2
My GDAL version is 2.2.4