texts:
I'm trying to write a simple qgis c++ application (Linux Mint).

My goal is to 'talk' to QGis programmatically in C++. Plugins seem at first over-complicated, as all I need for now is to extract geometry info (no need for GUI).

Below the code I'm stuck with - from Python examples it seemed to be a good first shot. It simply outputs 'oh'. I'm getting the same result (!layer.isValid()) running the code from this pdf, after small edits to update to qgis 2.0. Tried two .shp files - one of them in fact a bunch of files exported from qgis (.shp, .shx, .dbf..).

What would be possible causes for a non-valid layer?

int main(int argc, char ** argv) {

QgsApplication app(argc, argv, true);

puts ("im alive");

QgsVectorLayer layer ("/home/koscianski/Code/qgis/shapefiles/builtups.shp", "a name I like", "ogr");

//QgsVectorLayer *layer = app.addVectorLayer ("/home/koscianski/Code/qgis/shapefiles/builtups.shp", "a name I like", "OGR");

if (layer.isValid())
   printf ("Found %d", (int) layer.featureCount());
else
   puts ("oh");

return 0;
}


texts: How to create quadralinear polygons around a set of points

texts:
I have a point shapefile of locations of trees in an orchard that are almost, but not quite, on a grid. The orchard is divided into plots made up of groups of trees that either in a rectangle of rows and columns or as a series of 5 trees in a single row.

I would like to create a polygon coverage of the plots where the edges of the plot are halfway between each row of trees and halfway between the last tree of one plot and the first tree of the next, so that the entire area is filled with no gaps. At the edges, the plots should extend 1/2 the row spacing and 1/2 the tree spacing beyond the trees within the plot. Essentially a body-centered lattice, like:

Not all plots have the same number of trees along the rows. The point data are just the location of trees. I need to create the polygons and then attribute them with the data for each plot. Each set of rows has the same number of trees and number of rows for each plot. (note the rows are defined North to South by the way the trees were planted so you could think of them as columns).

Is there a way to create these polygons without individually measuring the locations of the corner vertices? Alternatively I would settle for rectangles but that wouldn't be as rigorous.


texts: Initializing Google Maps Streetview object looking at the position

texts:
I have a streetview included in my web site using google API to load it.

var theloc = new google.maps.LatLng(lat, lon);
var panoramaOptions = {
  addressControl: false,
  position: theloc,
  pov: {
    heading: 0, /* North */
    pitch: -40, /* Look down at angle 40deg */
    zoom: 0
  },
  visible: true
};
var streetview = new google.maps.StreetViewPanorama(document.getElementById("streetview"),
                     panoramaOptions);

This works, but I have one issue:

I am interested in loading the view by looking at the position I refer to, but when loading this it always points to north, and the location I specify is not necessarily in view.

If pointing to maps.google.com using URL to open a street view direclty in Google, I can simply omit the heading part of the panorama options parameter (cbp) like this:

'                                         \/ The heading value space                        '
https://maps.google.com/maps?layer=c&cbp=0,,,,40&cbll=58.9582956492561,5.62897047851361

This causes Google to automatically center the view on showing the position specified in cbll, and I don't need to know what heading to use.

Is it possible to do the same in the javascript API? If I try to omit the heading in the pov object, or initializing heading with null, the street view fails to load. I have also tried to set heading=-1, but that just points always to north. I also tried to not set the pov object at all, but that also always points to north (plus I don't get to set vertical angle 40 degrees down like I want).


texts:
I am trying to overlay features from a WFS using OpenLayers, my geometries in database are in EPSG:27700 and my map's projection is EPSG:3857. My mapfile returns data in EPSG:27700.

Is there any way I can transform points into the map projection?

Here is my code which is not working.

var wfs = new OpenLayers.Layer.Vector(
            'layername',
            {
                strategies : [new OpenLayers.Strategy.BBOX()],
                protocol : new OpenLayers.Protocol.WFS({
                    url : 'proxy.php?url=xxx',
                    featureType : 'typename',
                    "geometryName": "msGeometry",
                    "featurePrefix": "ms",
                    "srsName": 'EPSG:3857',
                    "maxFeatures": 1000,
                    "version": "1.0.0"
                })
            }
    );


texts: Selecting lines from first table based on bounding Box of second table in PostGIS

texts:
In PostGIS, I am trying to get the elements from first table(A) that intersect with bounding box of second table(B) . I get the error using the query below.
What would be the correct syntax.

SELECT A.* from A,B where ST_Intersects(A.geom,SELECT ST_Extent(geom) as table_extent FROM B);

  ERROR:  syntax error at or near "SELECT" LINE 1: ...from rivers,india
  where ST_Intersects(rivers.geom,SELECT ST_...
                                                               ^
  ********** Error **********

  ERROR: syntax error at or near "SELECT" SQL state: 42601 Character: 67


texts: Coded value domains are not showing in pop-up in JavaScript API webapp

texts:
I have a web app that was built with JavaScript API 4.10. One of my layers is a map image layer that includes many sub-layers with some fields that have domains. The issue is that in the pop-up only the domain codes show, but not the values. If I add the map service in AGOL or Portal the domains values show correctly.
Here is the code where this specific map image layer was added:
var Utilities = new MapImageLayer({
url: "https://domain.../arcgis/rest/services/.../Utilities/MapServer",
title: "Utilities - CADD Master Files"
});
Utilities.watch("loaded", function() {
    for (i=0;i<Utilities.allSublayers.length;i++) {
        var this_layer = Utilities.allSublayers.items[i];
        //this_layer.visible = false;
        if (this_layer.sublayers == null) {
          crtPopUpTemplate(Utilities.url, this_layer.id,this_layer);
        }
        else {
            //console.log("This layer has a sublayer");
            //In production mode, remove console.logs
        }
    }
  })

Code for the pop-up template function (used for almost all layers in the app):
function crtPopUpTemplate(basicURL, layerID, layerVName) {
    var url = basicURL+'/'+layerID.toString()+'?f=json';
    esriRequest(url, {responseType: "json"}).
      then(function(response){
        Alldata = response.data;
        var titleLayer = Alldata.name;
        if (titleLayer.includes('_') == true) {
            resultIf = titleLayer.replace(/_/g," ");
        } else {
          resultIf = titleLayer;
        }

        let fieldInfosArray = [];
        for(let i=0; i<Alldata.fields.length; i++) {
          fieldInfosArray[i] = {fieldName: Alldata.fields[i].name, label: Alldata.fields[i].alias}
        }

        template1 = {
        title: resultIf,
        content: [{
         type: "fields",
         fieldInfos: fieldInfosArray
        }]
        };
        layerVName.popupTemplate=template1;
      });
  }

Any ideas?


texts:
Using QGIS 3.4

I have a set of points and I am using the atlas tool to generate reports for each point.

I'm looking for a way to only create reports for points that intersect a polygon.

These polygons change all the time, so im looking for a more on the fly way to filter these points, rather than the traditional select and save as a new layer.

In the atlas tools I see there is a "filter with" option and it gives you the ability to write an expression.

How can I create an expression/function to do the filter.


texts: Is there a way to run spatial join within tolerance?

texts:
I have some data that should be spatially joined, but from small precision differences (less than 1 in) they aren't joined. I've found no way to add tolerance to qgis spatial join operation.
Is there a way to config a tolerance value, or a plugin to correct small distances?


texts: Cannot store raster in custom database schema

texts:
I want to store an elevation raster in a postgis database. However a problem occurs. I use the following command:

raster2pgsql -I -C -e -Y -F -s 4326 -t auto dem.tif schemaname.filename | psql -U blabla -d blabla -h localhost -p 5432

As long as I put the raster in the public schema, the command works properly. But when I change the schema name to some other database schema I have created, the command just returns ERROR: schema "schemaname" does not exist.

Do I overlook something?

Edit: the full error looks like this:

ERROR:  schema "schemaname" does not exist
ERROR:  schema "schemaname" does not exist
invalid command \.
ERROR:  syntax error at or near "0100000100"
LINE 1: 0100000100B5EDD83C22BEF33E0CC0E6158FFAF2BE5981483035A34FC08E...
        ^
ERROR:  schema "schemaname" does not exist
ERROR:  The value provided for schema is invalid
CONTEXT:  PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean) line 53 at RETURN


texts: I am using a Tile Layer for my map and XYZ for the source of this tile layer in OpenLayers. I am providing PNG files as Z/X/Y. I want to replicate the real world coordinates on my map. The files that i am providing is of a certain country. As they are PNG files they do not have coordinates like open source maps but I want to map the real world and my source. Is it possible to do this in OpenLayers?


texts: Delaunay triangulation algorithm in shapely producing erratic result

texts:
I would like to triangulate the following polygon:

The wkt for the geometry is given as follows:

POLYGON ((3 0, 2 0, 2 0.75, 2.5 0.75, 2.5 0.6, 2.25 0.6, 2.25 0.2, 3 0.2, 3 0))

I ran this code to generate the triangles

triangles = shapely.ops.triangulate(res_intersection.geometry.values[0])
triangles_gdf = geopandas.GeoDataFrame()
triangles_gdf.geometry = triangles

When I take a look at the resulting polygon this is what I get:

However this is not what I want. The algorithm seems to generate some extra triangles which Are not contained in the original geometry. Why does this occur? Is there a workaround for this? Ultimately I would like to triangulate a large irregular polygon with holes.


texts: Difference between two atmospherically corrected Images

texts:
I corrected Landsat 5 image with two different methods, first atmospheric correction through calibration (means converting raw image into radiance), second atmospheric correction through QUAC atmospheric correction. The results I got are little bit change, the first one is a little brighter than the second one, so what is the difference between these two methods? and why my Image outputs are like that.
Radiance Image

Output of QUACT atmospheric correction


texts: Does Configurable Map Viewer support layer transparency?

texts: Configurable Map Viewer 2.0.0 Beta 2 is not honoring the layer transparency of individual layers in the map service (as configured in ArcMap).  On the other hand, transparency for entire map services as configured in CMV configuration files works just fine, but insufficient.  Is this a bug, a version incompatibility (I am currently being forced to use an ArcGIS 10.3.1 server), or is there something I missed about configuring CMV?


texts:
I watched the tutorial "Create your own data using the CartoDB Editor" about adding a feature but in the editor, I couldnt seem to find the adding a line feature.

When I click the button of add feature, there were no options of add line, the default response was just to add a point


texts: How to correctly convert data from OSM with osm2po to use it with A* algorithm?

texts: I want to make my A-Star fast, by using heuristics. At the moment it is very slow and even slower than my Dijkstra. I thougt A-Star is in the worst case as fast as Dijkstra. I convert my OSM bounding box wiht osm2po to a postgresql database with postgis support and then I do my routing with java on android. Do I have to consider something in the osm2po configuration file to get the right weight? For the heuristics calculation I use the Manhattan distance.


texts: Does anybody know if Qfield can run on Ipad? My understanding is it was build to run on Android 4.3 or later version. Anybodu knows if it can run on with IOS?


texts: Export all coded value domains from a geodatabase

texts:
On the ESRI-L mailing list this morning there was a question about how to see or export all the coded value domains for a geodatabase. The goal is to present the content of the domains in a tabular form, so they're easy to read.

The DomainToTable tool does this easily for a single domain, but when there are many domains it quickly grows tiresome. The best advice I could give was to the batch processing feature, but even that requires knowing or looking up the names of domains individually.

Surely there's a better way?


texts: which line feature connected to which two point feature?

texts:
i have point Feature (man hole) and line feature (cable).
just i want to know which cable length is connected in between which two manhole. if i put spatial joint it will take only starting point or end point. i want know the both point.
Please revert me further clarification.
i am new to the GIS Field. so please kindly help me out.


texts: Deploying configurable application template for ArcGIS Online or Portal for ArcGIS?

texts:
I am trying to use a configurable application template, Interactive Filter (Live Example Here) to create an online map in ArcGIS Online to be shared with others. I want my users to be able to filter the different layers that I have on my map. Those layers are uploaded as hosted layers under "Contents" in my ArcGIS Online account.

The Interactive Filter configurable application template is on GitHub and I have download and un-zipped it.

I really don't understand where to begin because I am very new to interactive and online maps. As far as I understand from the "Readme" file, I have to configure the JavaScript files and HTML and then upload the app in my ArcGIS Online account and share the URL of the app to the end users.

I want my map to be Sweden not USA, and I have multiple layers on the Sweden polygon that I want to use the filter on.

Where do I start?


texts: Get SpatiaLite Database path from Vector Layer in PyQGIS

texts:
What if I would like to have the full SpatiaLite database path from a vector layer loaded in the TOC with PyQGIS?

With a very simple example, by storing active vector layer in a variable:

vl = iface.activeLayer()

one can retrieve the providerType and the dataSourceUri:

vl.providerType()
u'spatialite'

vl.dataProvider().dataSourceUri()
vl.dataProvider().dataSourceUri()
u'dbname=\'/home/matteo/myDataBase.sqlite\' table="river_layer" (geometry) sql='

Is there a simple method to get as a string, just the path of the database (so. /home/matteo/myDataBase.sqlite)?


texts:
I have a 25 large geoTiff files, representing a DEM, to be published via Geoserver in order to have OGC WCS.
The size of each uncompressed file ranges between 235 MB and 6GB.
What would it be the best method to follow considering performance?
I read a lot of documentation about image mosaic plugin and imagePyramid plugin, but I do not still figure out what to do, whether adding inner tiles and overviews on all geotiffs then publish using image mosaic or merging those file into one file then using pyramid plugin.


texts: Using ArcMap 10.3 with PostgreSQL 9.3.x database and PostGIS 2.1.1?

texts:
I have a PostGIS database that we use for data storage. I would like to open this data in ArcMap instead of storing it in another format.

In previous versions of ArcMap it has been possible to open data from PostGIS with varying combinations of ArcGIS and PostgreSQL as supported by ESRI.

I now have a PostgreSQL 9.3 database (v9.3.3, also another one which is v9.3.4) with PostGIS v2.1.1 installed.

In previous versions of ArcGIS I have been able to view PostGIS data. I haven't enabled ESRI geodatabse or anything on my database, it is a normal PostGIS database.

I upgraded to ArcGIS 10.3 and ArcMap will not open PostGIS layers.

ESRI state that with ArcGIS 10.3 the required files for using PostreSQL/PostGIS are included as standard.

I have checked that the dll files included with ArcGIS are indeed 32 bit in C:\Program Files (x86)\ArcGIS\Desktop10.3 (since ArcMap, ArcCatalog etc are 32 bit you need the 32 bit files) and they are. ArcGIS Pro installs with the 64 bit version in C:\Program Files\ArcGIS Pro.

I cannot find very much information on this hence the question here.

Anyone had any success with ArcGIS 10.3 (ArcMap not Pro) and opening PostGIS data? I am going to upgrade my database to PostgreSQL 9.3.5 since ESRI state that as the supported version (ESRI website), to check this isn't the reason for the error. Other than that I am out of ideas seeing as I do not get a helpful error message. The error message is as below:


texts: How to export Spatial data from MsSQL to MySQL using OpenSource/Free tools

texts: I have a few tables with geometry fields in MsSQL Server 2008. I would like to export the whole table with data to MsSQL Server. Are there any open source/free tools that does this job?


texts: Projection problem for Chinese river file into WGS84 43N

texts:
I am using QGIS 3.10 and having issues with (re)projection of two shapefiles. For a map I am trying to make of a river network in Xinjiang China and South East Kazakhstan, I have found two shapefiles of Chinese rivers.
These were originally uploaded to the site in "Albers Conic Equal Area, datum WGS84". My project is in WGS84 43N. Adding the layers initially gave the error of "invalid CRS".
Already tried reprojection, saving in different CRS, adding the shapefiles as a first layer in a new project - all to no avail. The shapefiles are now in the desired projection, yet they still hover above the South Indian Ocean.
Any tips on what I could try and do?


texts: Connecting Add Feature Class to Topology to Add Rule to Topology in ModelBuilder?

texts:
I have a question regarding to how to connect the outputs of "Add Feature Class to Topology" to "Add Rule to Topology" in ModelBuilder.
As the picture shows, I want add a rule called Must Cover Each Other (Area-Area) in this Topology, I did "Add Feature Class to Topology" for both feature classes before adding this rule, however, it just allow me input one of the 2 outputs as "Input Topology".
How do I connect the 2 outputs of "Add Feature Class to Topology" to "Add Rule to Topology" in ModelBuilder correctly?
P.S.: I also tried to connect the feature classes and topology directly to "Add Rule to Topology", it returns an error and failed...


texts:
I have a file with many polylines that I decoded with polyline decode on python. The result was a csv with this format:

[(-23.56632, -46.68435),
  (-23.56629, -46.68428),
  (-23.56589, -46.68346),
(-23.56307, -46.6786),
  (-23.56308, -46.67861)],
 [(-23.54487, -46.66172),
  (-23.54485, -46.66192),
  (-23.54464, -46.66237),
  (-23.54461, -46.66244)], ...]

Now I want to open this in QGIS to see the many routes I have but I don't know how to do this.

I've tried "Add delimeter text" > WKT  but it still shows the same error:

Layer is invalid: layer in file:///C:/Users/... is not a valid layer and can't be add to map

Then I tried to change my lat lon appear like this

link

but I still receive the same error. What can I do?


texts: Sourcing GIS data from First World War (WW1)?

texts:
Simple question, but does anyone know where I can get some detailed (eg 1:10m scale) vector data from the First World War?

I've been through the archives here, but it seems the best sites for it have already closed.


texts: Function to make a trajectory more straight

texts:
I'm new at PostGIS and I have to do a project for my university course.
I'm trying to analyze a trajectory of points that I downloaded from my runtastic app but it's a mess.
There are some points where the trajectory goes back and forth (green line in the image), while my real trajectory was straight (pink).

Is there a query or function I can use to delete those points and make my trajectory a bit more straight?


texts: Empty reply from server when downloading NLCD 2016 in R

texts:
I am trying to download the NLCD 2016 Land Cover Dataset in RStudio. The area I am trying to get data for is roughly that of New York State. The code I am using at this point is:
library(FedData)
library(sf)

NYSbb <- data.frame(lat = c(40,45.2),
                 lon = c(-71,-80))
NYSbb_df <- st_as_sf(NYSbb,
                    coords = c("lon","lat"),
                    crs = 4326)
NLCD <- get_nlcd(template = NYSbb_df,
         label = "NYS",
         dataset = "Land_Cover",
         year = 2016,
         force.redo = T) # Left this line in here because of cancelled download attempts

The error I receive after running this code is the following:
Error in curl::curl_fetch_disk(url, x$path, handle = handle) :
  Empty reply from server

I understand that a curl error is some sort of a server error. I am on a personal computer, and I have no problem accessing the NLCD website from my internet browser.
EDIT: Had label and dataset arguments switched, error persists.


texts: I am using ArcGIS 10.1 and I would like to change the coordinate system of a layer that is already attached to a drawing so that it matches the other layers.


texts: osm2pgsql cannot connect to postgis database, running on a Bitnami project installer

texts:
Based on this guide I am trying to import an .osm.pbf file into my PostGIS database, which is running on my Postgres Bitnami installation here:

[bitnami@ip-xxx]/usr/share/proj$ which postgres
/opt/bitnami/postgresql/bin/postgres

but I am getting an error, where osm2pgsql seems unable to connect to my database:

[bitnami@xxx]/home/.../geoserver_data/OSM$ osm2pgsql -E 900913 -d geodb -U geouser -W -S /usr/share/osm2pgsql/default.style washington-latest.osm.pbf
osm2pgsql SVN version 0.82.0 (64bit id space)

Password:
Error: Connection to database failed: could not connect to server: No such file or directory
        Is the server running locally and accepting
        connections on Unix domain socket "/var/run/postgresql/.s.PGSQL.5432"?

Yet, my Postgres install is listening to port 5432 under /opt/bitnami/postgresql/.s.PGSQL.5432 as follows:

[bitnami@ip-xxx]/home/.../postgresql/share$ sudo netstat -nlp | grep 5432
tcp        0      0 127.0.0.1:5432          0.0.0.0:*               LISTEN      13969/postgres
unix  2      [ ACC ]     STREAM     LISTENING     1876491  13969/postgres      /opt/bitnami/postgresql/.s.PGSQL.5432

It seems that osm2pgsql is having trouble finding my database. I know that I have my postgres server running, as shown here:

[bitnami@ip-xxx]/usr/share/proj$ ps aux | grep postgres
postgres   982  0.0  0.0 247712   468 ?        S    Oct21   0:17 /usr/lib/postgresql/9.3/bin/postgres -D /var/lib/postgresql/9.3/main -c config_file=/etc/postgresql/9.3/main/postgresql.conf
postgres   986  0.0  0.6 247856  7028 ?        Ss   Oct21   0:00 postgres: checkpointer process
postgres   987  0.0  0.0 247712   140 ?        Ss   Oct21   0:14 postgres: writer process
postgres   988  0.0  0.0 247712   104 ?        Ss   Oct21   0:12 postgres: wal writer process
postgres   989  0.0  0.1 248568  1088 ?        Ss   Oct21   0:12 postgres: autovacuum launcher process
postgres   990  0.0  0.0 103520   536 ?        Ss   Oct21   0:22 postgres: stats collector process
postgres 11111  0.0  1.1 171996 11876 ?        S    04:28   0:00 /opt/bitnami/postgresql/bin/postgres -D /opt/bitnami/postgresql/data
postgres 11113  0.0  0.0 171996   876 ?        Ss   04:28   0:00 postgres: checkpointer process
postgres 11114  0.0  0.1 171996  1696 ?        Ss   04:28   0:00 postgres: writer process
postgres 11115  0.0  0.0 171996   932 ?        Ss   04:28   0:00 postgres: wal writer process
postgres 11116  0.0  0.2 172844  2084 ?        Ss   04:28   0:00 postgres: autovacuum launcher process
postgres 11117  0.0  0.1  27812  1120 ?        Ss   04:28   0:00 postgres: stats collector process
bitnami  13358  0.0  0.0  11752   924 pts/6    S+   09:10   0:00 grep --color=auto postgres
root     26079  0.0  0.0  64536    32 pts/7    S    Nov26   0:00 su postgres
postgres 26080  0.0  0.0  21100   200 pts/7    S    Nov26   0:00 bash

I also have the following under my /etc/postgresql/9.3/main/pg_hba.conf, /opt/bitnami/postgres/share/pg_hba.conf and my /opt/bitnami/postgres/data/pg_hba.conf, which seems like is should allow access for my geouser:

# TYPE  DATABASE        USER            ADDRESS                 METHOD
local   geodb           geouser                                 md5

Any suggestions about how I can import my OSM data?


texts:
I am using the current Qgis 2.6.1 and having problems using the measuring distance function. The 'Measure Line' is giving me distances which is completely wrong and the scale shown was too big also. The message displayed was that the map "...defined by Geographic Coordinate System but the map extent suggest it is actually Projected Coordinate System...". I tried rectifying it using the Project Properties menu but still no good.

It was alright when I used  the same data in the previous version (2.6.0)
The CRS = WGS 84.


texts: Assigning non-geographic data to table in MapInfo?

texts:
I'm trying to update a non-georeferenced table (Table A) to make it georeferenced/mappable.  I'm trying to do this using a common set of data also found in a georeferenced table (Table B).

So both tables have one column which contains the same information.  This information is a list of sites/polygons.

Table B has one row for each site while Table A has many rows for each site.

I need to make Table A's list of sites georeferenced or add a georeferenced column.  For example, if row 1 in Table A contains the site "X" ID then assign that row to polygon "X" in the Map Window.


texts: How to automatically get/assign contour height attributes?

texts:
I have quite a simple problem yet I can already guess the solution is going to be pretty complicated.
I'm using QGIS and was wondering if there is a simple way for the computer to automatically assign height attributes to contours. The .kml file does not include height data and I don't have a .DEM file, but my source does have a .pdf to give a general idea of the contour layout. You can find it here: http://www.carpworld.co.za/index.php/dam-contours/boskop-dam. Say I set the height of the outside(highest) contour and the most inner(deepest) contour of the lake as well as the contour interval, I can see the plausibility.
This is how the lake and the Attribute Table looks:

The features displayed in  Attribute Table are points and the contour lines are not joined yet (One contour is not displayed as as one feature) they are still segmented lines and not polylines, so this is something to keep in mind.

So to summarize: I need a simple quick way to join the contours and assign height values to them in a somewhat automated process.


texts: Plot in rstudio a NetCDF File - latbox by lonbox

texts:
I have two NetCDF-Files. One of them holds information about the temperature in 2m (t2m) , the other one about outgoing long wave radiation (olr). I want to make a scatter plot. The first part of this exercise is already done: the task was to select a lon-lat-box from each NetCDF-file. It's going from 50N to 70N with -180W to 180 E. That's where the "50N70N comes from in the file name). Both files have a lot of values.

It's the first time I am working with the package ncdf in Rstudio. So what I already have done was to get all the variables from the files. As there is only one variable within the files, this was no problem at all.

Now comes my question: I just tried to plot the two files as a scatter plot, but I have the feeling that the program is plotting for each point all temperature and radiation values. I have too many points. Of course this is not what I want. I want that every box has only one t2m to one olr value.

Can you guys help me how to say in the script that it should only plot one temperature value to one radiation value? How can I say to plot separately box by box ? I hope you guys understand my question. I will upload the current plot.

Maybe I am all wrong and the program is already doing it like I want.
Here are the print.ncdf from both files:

"file olr_50N70N.nc has 3 dimensions:"
"lon   Size: 144"
"lat   Size: 9"
"time   Size: 366"
"------------------------"
"file olr_50N70N.nc has 1 variables:"
"float var179[lon,lat,time]  Longname:var179 Missval:1e+30"

"file t2m_50N70N.nc has 3 dimensions:"
 "lon   Size: 144"
 "lat   Size: 9"
 "time   Size: 366"
 "------------------------"
 "file t2m_50N70N.nc has 1 variables:"
 "float var167[lon,lat,time]  Longname:var167 Missval:1e+30"

Here is my script:

library(ncdf)

library(maps)

t2m_50N70N <- open.ncdf("t2m_50N70N.nc")
float_t2m_50N70N <- get.var.ncdf(t2m_50N70N,"var167")
print.ncdf(t2m_50N70N)

olr_50N70N <- open.ncdf("olr_50N70N.nc")
float_olr_50N70N <- get.var.ncdf(olr_50N70N,"var179")
print.ncdf(olr_50N70N)

plot(float_t2m_50N70N, float_olr_50N70N, type="p", pch=20, col="red", lty=3)


texts:
I would like to use Osmfilter for extracting specific relations by ID (and their dependent ways and nodes).

All I've found was filtering by tag/key.

Is there a way to filter by ID?


texts: About event to Latitude and Longitude in C#

texts: I am beginning to study Ready GXP and I would like to know which event receives Latitude and Longitude and passes them to the browser for execution in C#.


texts: Converting AutoCAD mtext & mleaders into multiple GIS attributes

texts:
I have a large utility system in AutoCAD dwg format.  I need to convert it into ArcGIS shp or gdb format.  It consists of structures and pipes all with associated text (mtext & multileader objects) all containing several parameters that will need to each broken out into it's own attribute field (eg: one attribute field for material, one for size, one for top invert, etc etc, see link for example)

Is there any set of FME/AGDI transformers or combination of that can automate or partially automate this process by reading from the text and converting it to multiple GIS attributes?  I've found some transformers that can read & hack up the text, but I've found nothing that looks like it can associate the text with its respective features.  There are over ten thousand features that need attribution and my aim is to avoid having to look individually at each feature to populate its fields.  I have AutoCAD Map3d, ArcGIS and it's Data Interoperability extension at my disposal and can get the full FME program if needed.


texts: can I turn off zoom resolution in openlayers 3 drag zoom

texts: Using an openlayers 3 map. I want to do drag zoom to in-between resolutions rather than to the nearest integer zoom.  I can see that pinchzoom has a constrainResolution but there doesn't seem to be anything on dragzoom. Openlayers 2 seemed to have a fractionalzoom that did what I want, can I get the same effect in openlayers 3?


texts: Converting a 32-bit GeoTIFF to map tiles using gdal2tiles

texts:
I have a 32-bit GeoTIFF file I'm trying to convert to map tiles using gdal2tiles.py.

What the GeoTIFF is supposed to look like when loaded straight into QGIS (no scaling, styling, etc):

And here is the output of gdalinfo input.tif:

bash-3.2$ gdalinfo input.tif
Driver: GTiff/GeoTIFF
Files: input.tif
       input.tif.aux.xml
Size is 294, 192
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-83.400000000000006,35.200000000000003)
Pixel Size = (0.016666666666667,-0.016666666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -83.4000000,  35.2000000) ( 83d24' 0.00"W, 35d12' 0.00"N)
Lower Left  ( -83.4000000,  32.0000000) ( 83d24' 0.00"W, 32d 0' 0.00"N)
Upper Right ( -78.5000000,  35.2000000) ( 78d30' 0.00"W, 35d12' 0.00"N)
Lower Right ( -78.5000000,  32.0000000) ( 78d30' 0.00"W, 32d 0' 0.00"N)
Center      ( -80.9500000,  33.6000000) ( 80d57' 0.00"W, 33d36' 0.00"N)
Band 1 Block=294x3 Type=Float64, ColorInterp=Gray
  Min=0.005 Max=0.759
  Minimum=0.005, Maximum=0.759, Mean=0.057, StdDev=0.076
  NoData Value=-1.69999999999999994e+308
  Metadata:
    STATISTICS_MAXIMUM=0.75856912136078
    STATISTICS_MEAN=0.057134358222004
    STATISTICS_MINIMUM=0.0052262963727117
    STATISTICS_STDDEV=0.075826669687323

What a sample tile looks like when converted using gdal2tiles.py input.tif:

OK, so clearly there's something going on with floating points and possibly number signing. So I try to convert it:

gdal_translate input.tif output.tif -ot UInt16

... And then convert to tiles:

Result: Only values over 0.5 carry over as black. Everything else is white.

I feel like there is something obvious that I'm missing. Does anyone have any insights into what's going on here?

Things I've done:

-Try all the different convert options:

Byte: Monochrome conversion

Int16: All black

UInt16: Monochrome conversion

UInt32: Monochrome conversion

Int32: All black

Float32: TV snow noise

Float64: TV snow noise

CInt16: All black

CInt32: All black

CFloat32: TV snow noise

CFloat64: TV snow noise


texts: How to find the distance to the nearest coastline from a land-based latitude/longitude point?

texts:
I'm completely new to GIS and unsure what I need to solve this problem:

In the image above turquoise is actually land. Sea is white.

I tried using the NASA Distance to the Nearest Coast data dump by loading it with NumPy and performing a binary search, which worked and was blazingly fast, but the dataset's resolution wasn't high enough. I need a resolution of about 10 metres.

I've now downloaded the GSHHG shapefile (shown above loaded in QGIS), but I'm unsure what to do with it. I'm seeing a lot of posts here mentioning PostGIS and SpatiaLite, but I'm not sure if these would be fast enough, or even the right tool. The maximum lookup time needs to be under 100ms as I will be running this against 10,000s of locations during the back-end indexing phase of a travel search engine website.

EDIT:

I'm trying to generate the raster from the shapefile, but I'm not sure what resolution to use for the -tr option. According to https://gis.stackexchange.com/a/8655/30920 I need 0.0001, as 11 metres would be ideal.

Below shows the output of -tr set to 1, 0.1 and 0.01 respectively. The file size seems to be growing by 100x for each decimal place. With 0.01 being 4.1gb that would make 0.0001 41tb (way too big).

I'm running:

gdal_rasterize -a level -tr 0.01 0.01 -l GSHHS_f_L1 GSHHS_f_L1.shp 0.01.tif

Is this right?


texts: GEE: Batch Download S2 Image Collection with pixel values = NoData

texts:
I am downloading the S2 image collection with all the bands(Band 1-Band 12)based on multiple boundaries.  The code is working well and it manages to download the GeoTiff S2 image with all Band 1-Band 12.  However, all the downloaded images have pixel values of NoData when I opened them in ArcMap, QGIS and SNAP.  I try to calculate the statistic for raster but not work too.  What did I do wrong?
Here's my code https://code.earthengine.google.com/41f205709fd0cda8e313d888e47317fb


texts: How to reproject Raster/Vector in Global Mapper with custom geographic parameters?

texts:
I am using Global Mapper 21 and i want to reproject raster/vector data with custom geographic parameters.

I can do it with ArcGIS (create custom geographic transformation). I have 7 parameters that can be used on method coordinate frame (at least in this name i can find this method in ArcGIS) but cant find any option to use this parameters in Global Mapper. I want to transform from epsg:3127(or 6870) to 32634 and vice versa also from 28404(or 2462) to 32634, i have 7 parameters already.


texts: Display WFS (GeoJSON) from GeoServer with OpenLayers on Django

texts:
I want to display point data from PostGis to my Django application. I created WFS service from GeoServer and want to display using OpenLayers. I went through many posts but none of them worked
My code is
myproject/myapp/templates/map.html

    <link rel="stylesheet" href="http://openlayers.org/en/v3.12.1/css/ol.css" type="text/css">

Plant Functional Types

<button type="button" onclick="go()">Actual plant functional type</button>
<button type="button" >Potential plant functional type</button>
<button type="button" > Plant functional types likely to change</button>
<button type="button" > Upload species data</button>

when I run localhost:8080/pft
only map gets dispalyed. Features are not displayed. Please help!

geoserver log shows this

2016-01-06 17:03:02,960 INFO [wfs.json] - about to encode JSON
2016-01-06 17:03:02,964 INFO [geoserver.wfs] - Request: getServiceInfo
2016-01-06 17:03:02,966 INFO [geoserver.wfs] - Request: getFeature service = WFS version = 1.1.0 baseUrl = localhost:8080/geoserver query[0]: typeName[0] = {PFT}pft_pft maxFeatures = 50 outputFormat = application/json resultType = results


texts: Determining best neighbourhood search radius for Topographic Position Index (TPI) for predictive soil mapping?

texts:
I'm not sure how to determine the best input values for a search radius to create a TPI for the purpose of predictive soil mapping.

Would this depend on the spatial extent of the area I wish to map?


texts: update pgrouting from 2.0.0 to 2.2.0 due to malfunctioning

texts:
I have postgresql 9.1 and postgis 2.0 on ubuntu 12.04 currently installed, but when I try to execute the following command:

psql -d routing_test
ALTER EXTENSION pgrouting UPDATE TO "2.2.0";

output is the following:

ERROR:  extension "pgrouting" has no update path from version "2.0.0" to version "2.2.0"

Here is a screenshot from a bug in version 2.0.0, which I hope to fix with updating:


texts: Is it possible to show data when hovering over an item with qgis2web (OL3)?

texts: I'm using the qgis2web plugin and I want to show the records behind the different items on the map when I am hovering over them. Is that possible when you have layers devided into groups? And my output file needs to be Openlayers 3 type.


texts:
I've used GDAL_translate with "--config GDAL_PDF_LAYERS" option to convert US Topo GeoPDFs to GeoTIFF and it works great for creating collarless topos based on the new format.

However, I can't get it to remove the "Images" layer from the US Topo PDFs for North Dakota (Fargo area) no matter what I try, short of flattening the PDF and discarding spatial component, which defeats the purpose of GDAL conversion. The US Topos for MN in the adjacent area convert as expected. It's something to do with the ND topos.

Has anyone come across this and how can I get GDAL to remove the "Images" layer from the ND US Topos?


texts: How to snap a GPX-track to a road network?

texts:
I want to snap the GPX-track to the road network as the image below shows. I found out, that "Hausdorff Distance" is exactly what I need and that it is usable with PostGIS. But unfortunately I have no idea how to work with it and even I don't know where to start. I'm not very familiar with PostGIS.

How do I use ST_HausdorffDistance (http://postgis.net/docs/ST_HausdorffDistance.html)? Is there any tutorial or workshop?

Im using qgis. The road-network comes from OSM (osm2pgsql).


texts: PolarStereographic cannot set latitudeOfOrigin (standardParallel) in Geotools

texts:
We have geotiff file with metadata as below, and we are not able to create the proper PolarStereographic coordinate system.

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",10],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

During transformation units from PolarStereographic projection to EPSG3857_WEBMERCATOR projection, standardParallel variable from PolarStereographic class is used. And despite fact, that latitude_of_origin is equal 60
(see configuration above), standardParallel is set automatically to 90 degrees. That's why projection after transforming is moved to the top.
We are using version 17.2

Is there a way to set latitude_of_origin of PolarStereographic projection, other than Java Reflection API?


texts: Creating Buffer from field in Attribute table

texts:
I have a shapefile in ArcMap with about 30 polygons and want create an inside buffer, using a different value for the distance of each buffer on each polygon. I was wondering if there was a way to use a field in the  attribute table of that shapefile to created an inside buffer, and if there is a way, how would I specify the units for that field in the attribute table.

Thanks!


texts: Shapefile created from Whitebox GAT has no *. prj file?

texts:
I am new to Whitebox GAT.

I found out that the shapefile I created in it does not have projection file.

Did I do something wrong?


texts: QGIS: Updating labels based on delimited text layer

texts:
I have a csv file that I am joining to a layer, then creating labels based on the joined layer.  When creating the delimited text layer, I include the option to watch the layer.  At the time of creation, the labels are correct.

When the csv file is updated with new values in the columns, the labels are not updated.  I have even tried closing the QGIS project and reopening and they still aren't updated.

Short of deleting the delimited text layer and re-adding a new delimited text layer and rejoining to my layer and redefining the label format, how can I automatically update the labels based on an updated delimited text layer csv file?


texts: How do you deploy a custom widget in Web AppBuilder for ArcGIS Online or Portal for ArcGIS?

texts: I am trying to deploy a custom widget, (Save Session), in Web AppBuilder for ArcGIS Online. I want the users of my application to be able to save the layers they have checked-on and -off in their session for use in a future session, but I've never installed a custom widget before.  Where do I start?


texts: How to draw route from OSRM on right road using OpenLayers?

texts:
Asking to OSRM for a route, OSRM is returning the geometry:

"route_geometry":"ity_rArqlxIvTv|AoCrY"

After decode them, I'm getting 3 GPS coordinates:

43.529557,          -5.659434
43.529208999999994, -5.660934
43.529281,          -5.66136

But with just 3 points, I'm getting this route :(

How can I draw a right route on the road with openlayers? Like this:


texts: Snapping doesn't work when many enterprise layers load

texts:
I am using Arcmap 10.6, when I have only limited number of enterprise layers loaded, the snapping experience process is painful. My cursor turns to "Busy" shape (the blue circle) and can't snap to the feature. I have no idea why it happens. Any solution?

Notes:
Everything is hosted on a server through a Virtual Machine. The virtual machine includes all the mapping software for an enterprise level, along with SQL 2012.

Mapping Software:  Version 10.6

ArcGIS Server

Portal for ArcGIS

Two ArcGIS Web Adaptors

SQL Server 2012

Operating System

Windows Server 2012 R2

Memory

32 gigs

Internet Bandwidth

50/50 mbps up/down

Tested:

Rebuilt the entire map..
Checked the map's data source in ArcCatalog
Ran MXD Doctor
Restarted ArcGIS services in Services
Set the Default map as Default_old and generated another Default map
Full Compression (Includes Reconcile & Post w/ Rebuild Indexes)
Python ~ nightly routine w/ log (Successful)
Restricted memory limits in SQL
Checked the snapping tool bar
Checked the tables
Restarted the database


texts: running scripts of codeblocks in arcpy with syntax error

texts:
I am having a trouble in running the following scripts:

import arcpy
arcpy.env.workspace = r'C:\MyData\my_geodatabase.gdb'
max = 123
min = 0
mean_s = 12.3
expression1 = "test(!mydensity!)"
codeblock1 = """def test(value):
   if value <(%mean_s% + %min%)/2:     # if value <(12.3 + 0)/2:
      return "cat1"
   elif value < (%max% + %mean_s%)/2 and value >=(%min% + %mean_s%)/2:    #elif value < (123 + 12.3)/2 and value >=(12.3+ 0)/2:
      return "cat2"
   """
a ="a"
b = "b"
arcpy.CalculateField_management(a, b,expression1,"PYTHON", codeblock1)

this python script returns with error :
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 2)
Failed to execute (CalculateField).

but the python works if I did the change in the code indicated after the (# ...)in the code above(in other words,if I change those %mean%,%max%,%min% to the real numbers, it works)

I guess it is because the codeblock in python wont recognize the max min mean_s defined at the beginning or something else goes wrong?


texts: ArcMap 10.1 SP Spatial Join Problem (lines within polygon)

texts:
I have a line layer (faults) and a polygon layer (fishnet grid).  The fishnet has a risk of earthquake attributed to each polygon.  I want to take that risk value and give it to the fault lines which pass through that polygon.  I am attempting to use the spatial join tool to do this.

I have the target layer as the fault lines and join layer as the risk polygons.  I have selected the 'Within' operator and one-to-one join operation and set the risk values to be added as a field to the new layer with the merge rule set to sum (I want to add the risk of all the polygons that the line passes through).

The tool works for lines which are completely within a polygon and returns the risk value correctly but any lines which pass through 2 or more polygons are returned with a null value instead of the summed risk.  If I uncheck the keep all target features then it leaves them out, implying that the tool doesn't think the lines are 'within' any polygons.

Is this an error in Arcmap? and how can I fix this?  Or is that not how the 'within' operator works for lines within polygons?  My assumption is that the 'Completely_within' operator should return what I am getting from the 'within' operator (also this is how it's described in the tool help).

Green yellow and red lines have had the risk values correctly attributed but as you can see its only the lines which fall completely within a polygon.  The Purple liens are the ones which aren't being computed with the summed risk but are returned with null.

A screen cap of the tool before running (the warning is just the because the file name exists since i took it afterwards...)


texts:
I have a problem with the examples in the book "Python Geospatial Development". I can't write to the PostgreSQL database polygons from TM_WORLD_BORDERS (Chapter 7, page 202 in book). Code:

src_file = os.path.join('TM_WORLD_BORDERS-0.3.shp')

shapefile = osgeo.ogr.Open(src_file)
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField('NAME').decode('Latin-1')
    geometry = feature.GetGeometryRef()

    wkt = geometry.ExportToWkt()

    cursor.execute('INSERT INTO countries (name, outline) '+
                   'VALUES (%s, ST_PolygonFromText(%s, 4326))' +
                   '', (name.encode('utf8'), wkt))

    connection.commit()

Only 30 out of 246 polygons are written. But the geometry.IsValid() function shows that 242 country polygons are valid.


texts: Ogr2ogr transform shapefile into GeoJSON without the geometry

texts:
I'm trying to transform a set of shapefiles (49) that contains nearly 2500 LineStrings each and also different columns of information.

So my intention is to obtain only this additional information from the shapefile without all the geometries. For this I tried to transform the Shapefiles in GeoJSON (to later use inside some scripts with PHP) with ogr2ogr.

The command I tried is:

ogr2ogr -f GeoJSON  <geoJsonFile> <inputShape> -nlt NONE

So far this has worked perfectly in a local environment where I have ogr2ogr version 2.2. But in the server that I need to run the scripts (that has ogr2ogr version 1.9) it seems like the -nlt option is beeing ignored, as the resulting GeoJSON files still contains all the geometries.

I tried checking the documentation but it doesn't seem to be an issue with the version as the -nlt NONE option seems to be available just fine.

Is there any other reason why this may not be working as expected? Or maybe another parameter to achieve the same goal?

As a side note, I can't update the version that is available at that server, so updating to a newer version is out of my scope.


texts: How to deal with Roundabouts in PgRouting

texts:
Can anybody guide me to to deal with roundabouts in pgrouting?? I have some roundabouts in my road network and when i try to make route from one point to another, it gives me wrong direction where it finds the shortest path because length on road is my cost.

In the picture below, blue cross is my starting point and other side is mt destination.

I am using simple pgr_dijkstra function and i am not considering any direction. boolean parameter for direction is set false. And there is no information about direction in my data. I just want to know what should I do if I want correct routing on my roundabouts.

thanks in advance


texts: How to extract a small area from a big GTFS feed?

texts:
Is there a tool that supports extracting a small region from a bigger GTFS feed? Google's Transitfeed library contains a Merge tool but no Extract/Clip tool, as far as I can see. Other libraries I've found can only read and plot data but don't support manipulation.

Background info: The Transitfeed Feedvalidator shows a lot of errors for my GTFS file. I want to extract my region of interest to assess the local quality.

Recommendations on Twitter pointed me to FME which has both a GTFS reader and a writer. My simple workflow consists of a GTFS reader (bottom left), Shapefile reader (top left, for the area of interest polygon), a clipper and a GTFS writer for the output (only writing the stops and shapes files for now but eventually I'll want all files):

The issue so far is that the output shapes content (orange points) does not contain all the input points (green) which are within the yellow area of interest:

(Data source: https://transitfeeds.com/p/helsinki-regional-transport/735/latest/download)

Using  "Create Aggregates" = "No", the clipped shapes seem complete

but the stops definitely are not


texts: Where does MapInfo store its error and log files?

texts: I am having  errors opening a table with spatial co-ordinates in Microsoft SQL and my question in gis.stackexchange is here. I want to know where are the error and log files of Mapinfo stored. I searched over the Mapinfo installation but was unable to find them.


texts: Label expression help with loop field and multiple fields

texts:
I have 13 fields I need to Label all in one Label.  But I dont want to label values with zero.  Need to have Label before each field example (Field is (C_B) and need to read C&B, 2).  Here is what i have so far but not sure if this is correct.

def findlabel ([VF_CL], [C_B], [VF_C_R], [VL_C_R_Unp], [MH_F_C], [MH_F_C_Unp], [Reset_MH_F], [Reset_MH_1], [RtoG], [RtoG_Unpav], [Voids], [Joints],[Drop]):
    # Set up the Headers
    Label = "VF CL " + [VF_CL], "C&B " + [C_B], "VF C/R " + [VF_C_R], + "  (not sure if this is correct)
    # Loop through Fields
    flds = [[VF_CL], [C_B], [VF_C_R], [VL_C_R_Unp], [MH_F_C], [MH_F_C_Unp], [Reset_MH_F], [Reset_MH_1], [RtoG], [RtoG_Unpav], [Voids], [Joints],[Drop]]
    for fld in flds:
        if (fld)>0:
            label =
        else:
            return None
    return Findlabel


texts: Legend items to appear / disappear simultaneaously with Layer turned on/off in Layer Control

texts:
I have a number of layers that can be turned on and off with a Layer Control.
I’d like the description of each layer to show in the legend if it's turned on.
And disappear again as soon as the layer is turned off in the layer control.

I am new to js and Leaflet.


texts: Adding a new field and calculating area using Python

texts: Using Python, how do you add a new field to your shapefile and calculate the area in km2?


texts: UploadServiceDefinition_server to ArcGIS Server not working on IIS

texts:
I want to write a web using MVC that uploads a .sd file to my ArcGIS Server using Python (arcpy). In my controller I have:

public string upload_sd()
        {
            ProcessStartInfo start = new ProcessStartInfo();
            start.FileName = "C:\\Python27\\ArcGIS10.2\\python.exe";
            string add = Server.MapPath("~/Scripts/Python/uploadsd.py");
            String[] args = new string[] { add};
            Debug.WriteLine(args[0]);
            start.Arguments = string.Format("{0}", args[0]);
            start.UseShellExecute = false;
            start.RedirectStandardOutput = true;
            string result;
            using (Process process = Process.Start(start))
            {
                using (StreamReader reader = process.StandardOutput)
                {
                    result = reader.ReadToEnd();
                }
            }
            return result;
        }

in uploadsd.py I have:

import arcpy
try:
    arcpy.UploadServiceDefinition_server(r'RAIN17242.sd', r'arcgis on 10.0.0.201_6080 (publisher).ags')
except Exception as e: print(e)

When I publish my web to the IIS server, script does not stop and when I check my ArcGIS Server Manager, nothing is published.

When I run the script using the command prompt, it works fine.


texts: Finding intersection points of two polyline features using ArcObjects?

texts:
I'm currently working with street intersections trying to find point(or points) where 2 given streets intersect(!disjoin). Streets are all Polyline features, and I have to work with 2 features at given time. I think this way it can be determined if they intersect(haven't yet tested it),

    public List<IPoint> RoadIntersections(IPolyline poly1, IPolyline poly2)
    {
        IRelationalOperator ro = null;
        ro = poly1 as IRelationalOperator;

        ISegmentCollection sc = new PolylineClass() as ISegmentCollection;
        sc.AddSegment((ISegment)poly2);

        if (ro.Disjoint(sc as IGeometry) == false) //if they intersect
        {
            //TODO: find intersection points
        }
     ...

but what I would like to do is to find point (or points) where those two polylines intersect.


texts: How to load a NetCDF-file into QGIS using Python?

texts:
I have a CF-conform NetCDF-file which can be loaded into Quantum GIS (Version 1.8) by using the "Add Raster Layer"-button and subsequently choosing the file and the raster layers of interest. This works without problems. However, when trying to load the same file by following these steps from the python cookbook, the layer remains invalid. I tried it using the console:

>>> fileName = "/path/to/file/myfile.nc"
>>> fileInfo = QtCore.QFileInfo(fileName)
>>> baseName = fileInfo.baseName()
>>> rlayer = QgsRasterLayer(fileName, baseName)
>>> rlayer.isValid()
False

When trying to add the layer to the map registry, it doesn't show up:

>>> QgsMapLayerRegistry.instance().addMapLayer(rlayer)
<qgis.core.QgsRasterLayer object at 0x3715f80>

I tried the same procedure with a jpeg with world file and it works.

The cookbook mentions that GDAL might not support some formats. I don't know how to check this in QGIS, but since it works manually, I don't think that this is the problem.

Maybe, I have to specify the raster layer of interest somewhere in the code?

Any help is highly appreciated.


texts: How do you split a feature using an existing line in QGIS?

texts: Is there a way to split a feature (polygon) using an existing line (selected polyline) in QGIS?


texts:
I have a merged shapefile with line geometry. Every line has a Z value written in the database. I want to write this value from the database to the geometry as Z-value.

The basefiles are S-57 DEPCTN layers merged together using vectormerge in QGIS. The Z value is stated in the column [VALDCO]. I want this value in the geometry as Z value.


texts:
I want to edit a polygon as that I cut off a part of it. I want to freely make a line that define where the polygon should be seperated. I have tried to useedit vertices but it does not seem to be a good tool for the purpose.

I use ArcGis 10.1


texts:
Could anybody give me some advice on which one should I use? OpenLayers or Leaflet

I will implement OSM on a web-app for taxi management.
The requirements for the map part are:

Tracking of the drivers cars using a tablet through my web app.
Conversion of a given address to lat-long.
Showing a route between start and destination lat-long.
Completely free of charge


texts: Creating File Geodatabase from KML using QGIS?

texts:
Is there a way to create a FileGDB from a KML?

I have tried with QGIS 2.18.7 but cannot figure out a way.


texts:
I am doing my GIS course term project in which I am integrating 3 kinds of map with index overlay method to see how much heat stress level is on my walking trail. To do it I have the tree canopy shapefile which is polygons with 100 meter distance to the trail and I should convert this vector file to raster to continue working on index overlay method.

The problem is I don't know how to rasterize the tree canopy vector file when I want to consider only the polygons which distance is not more than 10 meter to the trail.

I tried it with "Polygon to Raster" but I don't know what value field should I choose. Moreover if I add a field in the  attribute table of my tree canopy how can I specify t calculate the distance between polygons to the path line?


texts: GeoWebCache Miss: request does not align to grid(s) 'EPSG27700'

texts:
I'm using GeoServer - GeoWebCache with one layer group set to British National Grid (EPSG:27700). All layers within the group are the same projection - EPSG:27700

I've created a gridset for EPSG27700 and set the bounds to match the extent of my layer group.

I'm using OpenLayers 3 to fetch my layer using the following code:

var floodLayer = new ol.layer.Tile({
            extent: [96092.79481278273, 4921.6053531033685, 692596.2537952941, 665725.907765107],

            source: new ol.source.TileWMS({
                url: 'http://GEOSERVER/geoserver/Flood/wms',
                params: { workspace: 'Flood', layers: 'FloodRisk', tiled: true },
                serverType: 'geoserver',
                projection: 'EPSG:27700',

               )
            })

The extents above match the extents for the layer and the grid set. Hence the "Origin" should be the same - right??

So why do all my tile requests say "request does not align". The gridset is set to 256X256, I'm requesting 256x256.

What am I missing here?

So after 1 week of trying to resolve this issue I have to say I am not much further forward. I've set the same extent on the layers, gridset and openlayers request. I've set the same set of resolutions that are used by the OS base layer on the cached gridset, and also set these resolutions on the Flood layer overlay. Everything is the same and should line up, but I'm still getting 95% cache misses with the occasional hit when zoomed in.


texts: Direct Connection from Postgres to ArcGIS results to two feature classes per table

texts:
I'm using ArcMap (10.3.1) to connect to Postgres (9.4.5), which has Postgis version 2.2.0 installed. Create/Enable Enterprise Geodatabase wasn't used in this database. It's just a connection to the database. The connection itself works fine, but I get duplicate feature classes for each Postgres table. Only difference I've found is that first one has ability for address geocoding and second one looks quite normal. Also the first one has function Create Feature Class, and the second one Create Layer. The Postgres tables were created with FME 2015.1.

When dragging the feature class to ArcMap View, they both contain geometries and exact same attributes.

I checked that this shouldn't be related to Id, since I'm using objectid attribute, and ArcGIS recognizes it as Object ID data type.

View in the catalog:

Postgres tables:


texts: How to Batch Export MXD to GIF using Arcpy

texts:
I am trying to batch export several MXD files into corresponding separate GIF files. I am using the script below but am not having much luck and am getting an error: 'Invalid MXD filename'

import arcpy, os
Workspace = "J:\\IrrigatedLands"
arcpy.env.workspace = Workspace
mxdList = arcpy.ListFiles("*.mxd")

for mxd in mxdList:
    gifname = "J:\\IrrigatedLands_temp.gif"
    arcpy.mapping.ExportToGIF(mxd, gifname)


texts:
Since version 2.12.x GeoServer has the ability to proxy a remote Web Map Tile Service (WMTS). Is it possible to publish the remote WMTS as a WMTS on my GeoServer? If yes, how can I realize it in Geoserver?

One use case for cascading the WMTS is to hide the api_key of the remote WMTS request URL.


texts:
I am constructing a geopackage from vector layers that have values for number of years, each year is an attribute (1965, 1995, 2025, 2055, 2085). I have styled the layers and added the styles to the geopackage. In the layer_styles file, each of my files has many styles attributed to each layer (Figure 1).

When I load the package however, only the default style loads and other styles need to be loaded from the styles interface rather than being listed along with the default style as previously (Figure 2).
Is there a way to import the geopackage with all styles being active and selectable from the right click menu or will I have to load them each time from the geopackage for each layer?
As a follow-up (if this is not possible), what is the purpose of the f_table_name attribute in layer_styles if not to facilitate loading?
Maybe this is a bug?

Things I have checked/tried:

The styles are in fact in the geopackage
The layer does have its name listed in f_table_name for multiple styles
I checked for guides, but all of the ones I found only mention the default layer style, which is not ideal if I want to have many styles for a same layer. I figure there must be a reason that you are able to save multiple styles...

Alternatively, is there a better practice method for this?
I want to avoid having to save multiple instances of the same layer within the geopackage just for styling.


texts: Building NetCDF driver in GDAL within GRASS GIS

texts:
I am trying to import a NetCDF into Grass GIS. I am using MacOS and GRASS 7.8.4.
The grass documentation states that: Both NetCDF import and export in GRASS GIS are basically driven by the GDAL library. Note that GDAL library needs to be compiled with NetCDF support to enable GRASS to read and write this format.
The above linked page leads me to the following page with instructions for building the NetCDF driver. On there, it suggests that I need to have something called HDF4 or HDF5. The HDF5 appears the be the preferred choice, so I downloaded it, installed it, and started the configuring.
./configure CC=gcc FC=gfortran F90=gfortran F77=gfortran CPPFLAGS=-DpgiFortran --with-pic --with-hdf5=/Users/myusername/Downloads/hdf/HDF5-1.12.0-Darwin/HDF_Group/HDF5/1.12.0

However, it runs into the following error:
configure: error: Can't find or link to the hdf5 library. Use --disable-netcdf-4, or see config.log for errors.

To be fair I cannot manually find a file called H5Fflush in my newly installed HDF5 directory, but I am not sure what to do about that.
I should mention that I grabbed the netcdf-c-4.7.4.zip NetCDF library here, which I hope was the correct version, but I had trouble figuring that out.
Note that I do not know if I have to configure any of this in a particular directory as I have never had to build anything in a software package quite so manually before. Thus this is all happening in my download directory.
Further down the road, I don't know how I will get this built library into GDAL within GRASS GIS. I found a suggestion on how to use homebrew to create a new version of GDAL with the correct component, but I do not know how to do so within the GRASS GIS program.


texts: OpenLayers and GeoExt are displaying my lines features as point features

texts: I stored data (road network: line features) on OpenGeo Geoserver compilation and when I tried to display these in openlayers or GeoExt, they are shown as point features.


texts: Geokettle Data Input from Oracle Database

texts:
i just installed GeoKettle and am trying to do some first steps in this.
Can anyone help me with the data input step to an Oracle DB?

I tried the OGR Input:

But i seem to be unable to find the right syntax for the connection.
Can anyone help me with the Data Source (maybe giving an exampl)?
And what exactly do i have to fill in at "Layer Name"? The table name in my database? (i integrated the table , A37_SIREN_UTM in my data source in my example, but putting it in as the layer name didnt work too.)

Another option would be the Table Input:

The ODBC connection geodaten is defined on my computer.
But here too i am unable to connect to the DB.


texts: Executing Fill in ArcPy gives ERROR 999999?

texts:
I am having trouble using Fill in arcpy as part of a tool. It works fine as a script, but when I create a tool, it fails. To localize the problem I modified the Help stand-alone script associated with Fill. I also tried creating a blank raster outFill before the Fill command, but with similar results.

# Name: Fill_Ex_02.py
# Description: Fills sinks in a surface raster.
# Requirements: Spatial Analyst Extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/sapyexamples/data"

inSurfaceRaster = arcpy.GetParameterAsText(0)

# Set local variables
zLimit = 3.28

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

# Execute FlowDirection
outFill = Fill(inSurfaceRaster, zLimit)

# Save the output
outFill.save("C:/sapyexamples/output/TestFilled.tif")
~~~~~~~~~~

The error message is:
Traceback (most recent call last):
  File "C:\Users\William\Desktop\Pit_test\Tools\Test_outfill_tbx.py", line 22, in <module>
    outFill = Fill(inSurfaceRaster, zLimit)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 2271, in Fill
    z_limit)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 2267, in Wrapper
    z_limit)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 999999: Error executing function.
("esriDataSourcesRaster.RasterWorkspace") Failed to create raster dataset
ERROR 010302: Unable to create the output raster: c:\windows\system32\t_t428.tif
ERROR 010067: Error in executing grid expression.
Failed to execute (Fill).

Failed to execute (TestOutfill).
Failed at Wed Sep 11 17:18:08 2019 (Elapsed Time: 0.49 seconds)


texts: QGIS 3 Show Red Vertex Crosses in Edit Mode

texts:
I don't know if this is a change in QGIS 3 or if I managed to turn it off somewhere, but in previous versions of QGIS when a layer was in edit mode there were red crosses on all the vertices.

I would like to have that again, I'm using 3.6. I imagine it is in settings somewhere but I can't find it!


texts: How do I make a public download page for my CartoDB table?

texts:
I have some CartoDB data that I'd like to make available to people as shapefile, geojson and csv.

I just came across http://opensandiego.cartodb.com/tables/5037/public/map which has a great map preview, with downloads in different format and a cool table view. I have a number of public CartoDB tables, but I can't find how to share my data in that nice page. Is there a button to create the url? Or some way I can figure out what that URL should be?


texts:
Is the QGIS maptip html code capable of handling iframes?

I have a point shapefile of river gauge station locations.  I would like to be able to embed the realtime gauge information from Gauge Map for each gauge into the map.  The gaugemap website allows you to embed the information into other websites through an iframe (see screengrab below).  However, copying the iframe text into the html section of the maptips for the layer in QGIS doesn't produce anything.

Am I doing it wrong?  Is it possible to pull this into QGIS through the maptips or might there be another way?


texts: Open source software for 3 D vector editing on top of stereo images

texts:
I am looking for an open source software that could be used to digitize 3D vector objects from a stereo pair of remotely sensed images. Photogrammetry softwares are too expensive for my use (giving an introduction course to 60 students, which means only a dozen of hours every year).

I know a few software that can process stereo images to generate ortho-images and DSM (e.g. MICMAC, OTB...). However, I would like to be able to display an epipolar pair in 3D (e.g. with anaglyph or other 3D devices) and draw some 3D objects directly in 3D. Does it exists (or are you aware of future plug ins to existing OS photogrammetry softwares)?


texts: How to reproject a set of data in Python to 4326?

texts:
I need to reproject a set of data in a numpy array with a shape of (121, 160) with

crs='+proj=merc +units=meters +a=6370000.0 +b=6370000.0 +lon_0=-76.0100021362 +lat_ts=20.3180007935 +nadgrids=@null'

To crs="epsg:4326" in the location

Upper Left  ( -96.5979538,   5.2124481) ( 96d35'52.63"W,  5d12'44.81"N)
Lower Left  ( -96.5979538,  34.0865555) ( 96d35'52.63"W, 34d 5'11.60"N)
Upper Right ( -55.4220505,   5.2124481) ( 55d25'19.38"W,  5d12'44.81"N)
Lower Right ( -55.4220505,  34.0865555) ( 55d25'19.38"W, 34d 5'11.60"N)
Center      ( -76.0100021,  19.6495018) ( 76d 0'36.01"W, 19d38'58.21"N)

I am currently doing

    crs4326 = CRS.from_string('epsg:4326')
    affine = from_origin(x.min(), y.max(), xLength/xCellCount, yLength/yCellCount)

    dstTransform, width, heigth = calculate_default_transform(
        originCRS,
        crs4326,
        xCellCount,
        yCellCount,
        x.min(),
        y.min(),
        x.max(),
        y.max()
    )

    destination = np.zeros((heigth, width), np.float32)
    reproject(
        var,
        destination,
        src_transform=affine,
        dst_transform=dstTransform,
        src_crs=netcdfAttributes['crs'],
        dst_crs=crs4326,
        resampling=Resampling.nearest
    )

But for some reason the image deforms look at the bottom


texts: Transforming admin units into hexagons: the assignment problem

texts:
I have an administrative tessellation with about 8000 units (the Census tracts of California) and, for the sake of visualisation, I would like to transform it into regular hexagons (or squares). I found the nice R package geogrid that works well on small cases (e.g. up to 300 units). To match the geometries is a O^3 or O^4 problem: The Hungarian algorithm would take way too long on even moderately large input (several hours or even days for 5000 objects).

Is there any other package (in R or on other platforms) to do this at a larger scale? Is there a sub-optimal alternative to the Hungarian algorithm that would solve this assignment problem?

This is an example from geogrid that shows exactly what I want to do:


texts:
Here is a leaflet map in R

https://lawsblog.netlify.com/post/leaflet-map/

The white lines are neighbourhood boundaries.
The green polygons are parks.
The parks are drawn over the boundaries obscuring the borders.
I want the white neighbourhood boundaries to be drawn on top of the parks so that the entire border is seen.

This sounds like an issue of layer ordering.

I thought I could use zIndex with leaflet to solve the problem. I've tried hoods = 200, parks = 300 and hoods = 200, parks = 100. I can't get it to work.

The relevant code:

# Add neighbourhoods
  addPolygons(data = hood_shp,
              fillColor = ~pal(be_per_cap),
              weight = 2,
              opacity = 1,
              color = "white",
              dashArray = "0",
              fillOpacity = 1.0,
              options = list(zIndex = 200) %>%
# Add parks
  addPolygons(data = parks,
              color = "green",
              weight = 1.0,
              opacity = 1.0,
              options = list(zIndex = 300),
              fillOpacity = 1.0)


texts: Creating Metes and Bounds text from a polygon in ArcGIS

texts:
Lots of questions about turning metes and bounds into a polygon, but how about going in reverse? Are there any existing tools or do I need to look at building a script myself?

Ultimately I'd like to generate points with x/y and at all the vertices and line segments with distance and direction as attributes (so I can label them).

This for submission to our local government as part of their legal requirements.


texts: Detect one-way road from origin_id and destination_id

texts:
How do I detect the one-way roads?
I have a graph G( MultiDiGraph(G) with origin(Id) and destination(Id) nodes. I would like to detect one-way roads. Example if there is an edge from the origin "1" to destination "2" and another edge from the destination "2" to the origin "1", then this is not a one-way road. And a column should be added to the edges table of this graph called oneway (boolean).
 G = nx.from_pandas_edgelist(edges,'origin', 'destination',
 ['key','id_edge', 'name', 'lanes','length', 'speed', 'capacity',
 'function', 'oneway', 'geometry'],
                              create_using=nx.MultiDiGraph())

For detection one-way roads, I have tried this (but doesn't work):
oneway = []
for u, v, data in G.edges(keys=False, data=True):
    for a,b, data2 in G.edges(keys=False, data=True):
        if u==b and v==a:
            oneway.append(True)
        else:
            oneway.append(False)

For the column that should be added, I tried this:
edges['oneway']=edges.apply(lambda x: False, axis=1)


texts: Using QGIS, how can I add a bounding box attribute?

texts:
I have a shapefile where I am exporting the features to GeoJSON objects with ogr2ogr based on a simple where clause. One thing that I would like to be included in the GeoJSON object is the  'bbox' member (outlined here) , is it possible to add this as an attribute to the shapefile using QGIS?

I've done something similar with centroids and added the x/y as separate attributes. Are there any other ways to programmatically add this to my geojson files?


texts:
Has anyone else noticed a problem getting OpenLayers to work within jsFiddle?

I can't see what's wrong with http://jsfiddle.net/dPdhZ/

or the others I've mentioned here:

https://stackoverflow.com/questions/10075862/trouble-getting-openlayers-examples-with-jquery-working-in-jsfiddle

But it's frustrating when they work fine as a static HTML document but not in a social, sharable format like jsFiddle.

EDIT: Don't assume that  works and be cautious with using 100% widths/height on CSS for the map in jsFiddle.


texts: I have a spatiallite database. In this database I have some attributes. One of them is the attribute Nr_. In this attribute the value must be assigned automatically and is not allowed to be a double. I would like it to be named as 1001,1002,1003 etc.


texts: Output in GeoJSON data format using Ruby

texts:
I have some data with Ruby script (obtained by geolocating, loading from config files etc) that I want to output as a GeoJSON file.
I have some ugly code that generates valid geojson files - but is there existing standard for doing this?
I need it into another script and thought about extracting my code into a gem. But I suspect that there is existing library for that, and I simply failed to find it.


texts: Debugging Serious Application Error from Python Toolbox?

texts:
I have written a few custom Python Toolboxes that my office uses.

One of my coworkers has an interesting issue that I cannot replicate on a different machine:

When he runs the tool, ArcGIS Desktop crashes and tries to send an error report.  The script does not begin running.  This happens when running the script in both ArcCatalog and ArcMap.

Only one of the scripts I have written causes this, others behave normally.  The script began causing this behavior recently, it previously worked fine for him.

Is there anything I should check?

Re-installing ArcGIS Desktop seems like the only thing that might fix this problem.


texts: Table Of Content widjet in arcgis javascript api

texts:
I'm using TOC widjet in my application using arcgis javascript api, it is working, but the problem is that the plus and minus images which collapse or expand the tree is not appearing.

I noticed that it comes from css file called carlo.css, I included the file, but the elements didn't take properties from this file, I'm not sure what's the problem, anyone can help??


texts:
I download a graph and plotted an accessability polygon on the projected graph according to tutorial

G = ox.project_graph(G) #

Then I wanted to project my accessibilty polygon g1 back to 'epsg:3857' with this code

import pyproj
from shapely.ops import transform

project = pyproj.Transformer.from_proj(
    pyproj.Proj('epsg:4326'),
    pyproj.Proj('epsg:3857')) # destination coordinate system
g2 = transform(project.transform, g1)  # apply projection

But when I print g2 there emerges an error despite the fact that g1 can be plotted and gives no errors

ValueError                                Traceback (most recent call last)
~\AppData\Roaming\Python\Python38\site-packages\IPython\core\formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

~\Anaconda3\envs\dsenv\lib\site-packages\shapely\geometry\base.py in _repr_svg_(self)
    392             if xmin == xmax and ymin == ymax:
    393                 # This is a point; buffer using an arbitrary size
--> 394                 xmin, ymin, xmax, ymax = self.buffer(1).bounds
    395             else:
    396                 # Expand bounds by a fraction of the data ranges

ValueError: not enough values to unpack (expected 4, got 0)

<shapely.geometry.polygon.Polygon at 0x18625802700>

Is this a correct way to convert a polygon? Or have I mistaken the projection type?

UPD
If I exchange the sequence of projections pyproj.Proj('epsg:4326'), pyproj.Proj('epsg:3857') to pyproj.Proj('epsg:3857'), pyproj.Proj('epsg:4326')
then I get the polygon but its centroid contains wrong cordinates. So maybe the correct question is, which projection is used for ox.project_graph?


texts: How do you set an AddIn button's Help source and topicID values at runtime?

texts:
I would like a way to dynamically set the help source/topicID for a button.

Example Snippet:

  <Button
    id="Acme_FullExtent"
    class="FullExtent"
    caption="Full Extent"
    tip="Zooms the map the the full extent."
    category="Acme"
    image="Images\FullExtent.png">
    <Help source="http://www.Acme.com/help.htm" topicID="SectionB"/>
  </Button>

More background: I need to link the button's F1 help to online Help that is published by someone else and which is a configuration nightmare for various levels of testing (e.g.Test, UAT vs Production). I like the idea of loosely coupled online AddIn help content, but the help content still needs to be tightly coupled with the appropriate environment (which is defined via configuration at runtime).

I suppose the fallback plan is to compile the help into a chm (like pre-10.1 AddIns) and carry that onboard the AddIn.


texts: Multiple symbols in the same points vector layer (data defined marker)

texts:
I've categorized multiple points for schools that are to be developed, reconstructed or demolished, with green, yellow and red colors if the schools capacity (the number of children) is increasing, staying the same or decreasing. The size of each point is data defined depending on how large the change in capacity is.

However, for some of the schools the changing capacity is unknown, and therefore 'NULL' in my attribute table for these schools. I'm looking for a way to give these schools a Star symbol instead of the circles, so that their size is not to be confused with the schools where the change in capacity is known.

QGIS, v 2.14. In Layer properties - Style I've categorized them as stated earlier. To give the ones with a known change of capacity a circle and the ones with unknown change of capacity (= NULL) a star symbol i tried the following so far with no success: in Layer properties - Style - Symbol Change-button, I added 2 Simple markers. 1st one a circle and 2nd one a star. In each one I then went to Data defined override - Edit, command: "Change_cap" is NULL and in the second one: "Change_cap" is NOT NULL However, it's still showing both markers instead of just one of them (depending on if change of capacity is known)


texts: Installing double MapGuide in same computer?

texts:
Can I install two MapGuide Open Source with different version in the same computer?

For example, I have installed MapguideOpenSource 2.6 in my computer, then I want to install another version of Mapguide that is MapguideOpenSource 3.0.


texts: Calculating how many acres fall within each 1 square mile section?

texts:
I am using ArcMap. I have a polygon shapefile of parcels that receive surface water from a water district. I also have a polygon shapefile of section boundaries. I want to calculate how many acres fall within each section (1 square mile). I also want to set my symbology to represent how many acres that receive water are in each section. Should be 0-640.
I've tried "Near" and "Intersect"


texts: How to find the area of polygon in geodjango?

texts:
I need to calculate the area of the Building (which the coordinates represents). I used transform() to calculate the area, but the calculated area is not accurate. Please help me calculating the area using the lat/lon coordinates.

          Class Parcel(models.Model):
              geom = models.MultiPolygonField(srid=3492)
              street_nbr = models.CharField(max_length=80)

          Class Primary(models.Model):
              geom  = models.MultiPolygonField(srid=4326)

      >>>parcel = Parcel.objects.filter(street_nbr__contains="2730", street_nam__contains="11TH AVE",city__contains="SACRA").values()
      >>>parcel[0]['geom'].geojson

     >>>build = Primary.objects.filter(goem__contained= parcel[0]['geom'].geojson).values()

     >>>build[0]['geom'].geojson
     '{"coordinates": [[[[-121.476756, 38.542154], [-121.47686, 38.542152], [-121.47686, 38.542136], [-121.476904, 38.542135], [-121.476904, 38.542132], [-121.476962, 38.542131], [-121.476966, 38.542295], [-121.476893, 38.542296], [-121.476892, 38.542278], [-121.476845, 38.542279], [-121.476845, 38.542281], [-121.476759, 38.542282], [-121.476756, 38.542154]]]], "type": "MultiPolygon"}'

      >>>x = build[0]['geom']
      >>>y = x.transform(3492,Clone=True)
      >>>y.area
      1636.8046

Here i'm getting the area as 1636.8046, I was supposed to get 2,327 sqft . I'm not able to find what I am doing wrong.


texts:
I have the following code in Python for ArcMap to first select a property polygon and then with this selection to intersect the layer which contains lines (sewers) (which pass through or are within this polygon). In the Attribute table of the sections is a column (called "BG_HAusPfad") with the file path to its corresponding PDF.

The code then takes this file path and appends the PDF somewhere else.
(Because we have all PDFs in one folder but need to do a folder per property with the PDFs of the sections of that polygon).

The code works for one polygon. The issue is that with the second property polygon the internal memory is still the same from before and appends the PDFs from the first property.

How can I overwrite this internal memory?

Here's the Code: (R1_Nummer is the number of the property polygon, aka its name. The folders of the properties are called this way)

import arcpy, os
from arcpy import env
arcpy.env.overwriteOutput = True
R1 = r"R1_NUMMER"
field = r"BG_HAusPfad"
# that's the PDF name
PathPDF = r"I:\…\All-Separate-PDFs"
MainFolder = r"I:\…\All-Properties"

mxd = arcpy.mapping.MapDocument("CURRENT")

for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = pageNum
    pageName = mxd.dataDrivenPages.pageRow.getValue("R1_NUMMER")
    Subfolder = os.path.join(MainFolder, str(pageName))
    arcpy.SelectLayerByAttribute_management(r"Property-polygons", "NEW_SELECTION", "{} = '{}'".format(R1, str(pageName)))
    arcpy.MakeFeatureLayer_management(selection, "in_memory\PropertySelected_lyr")
    arcpy.Intersect_analysis([r"Sewer-lines", "in_memory\PropertySelected_lyr"], "in_memory\Sewer_lyr", "ALL", "", "LINE")
    Cursor = arcpy.SearchCursor(r"in_memory\Sewer_lyr")
    for row in cursor:
            Protokolname = row.getValue(field)
            ProtipdfPfad = os.path.join(PathPDF, str(Protokolname) + ".pdf")
            try:
                HprotipdfName = os.path.join(str(Subfolder), str(Protokolname) + ".pdf")
                HprotipdfDoc = arcpy.mapping.PDFDocumentCreate(str(HprotipdfName))
                HprotipdfDoc.appendPages(str(ProtipdfPfad))
                HprotipdfDoc.saveAndClose()
                del HprotipdfDoc
            except:
                continue

So in the picture you see the attribute table of the sewers (the layer is actually called "A1 aufgenommen" (since it's german and very specific i just changed it in the Code to sewers for you)


texts: Failed to create reader using geoserver Image Pyramid plug-in

texts:
I am following instructions on this link,
http://docs.geoserver.org/stable/en/user/tutorials/imagepyramid/imagepyramid.html

I am pointing the URL to the folder normally but I get the following error,
"Could not list layers for this store, an error occurred retrieving them: Failed to create reader from file:///mnt/MapData/bmpyramid and hints null"

The difference is that in my case, I used a mosaic.vrt file to create the image pyramids. If that was the problem, what would the solution be?


texts: Calling for layer in ArcPy with ArcGIS Pro

texts:
Every single example I see online, the code has the mapname.listLayers() in order to call the layer.

Don't we have a way to call the layer by its name instead? The reason I'm asking this is because this code does not work at all:

aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.listMaps("Map")[0]
aprxMap.addDataFromPath(outfcname)
install_layer = aprxMap.listLayers(installname)
layer_simbology = r"C:\ArcGIS\LAYER\symbologylayer.lyrx"
arcpy.ApplySymbologyFromLayer_management(install_layer, layer_simbology)
aprxMap.removeLayer(install_layer)
aprx.save()

The code does not update the symbol, neither remove the layer from the map. But add the Feature Class (outfcname) to map.


texts: Definition Query expression limits in ArcGIS 10.2.2

texts:
I have a large database in ArcGIS 10.2.2, with over 2 500 000 point vectors. Each point has a unique number code.
I want to spatially subset my data, so that I can extract only points with a specific codes.
I tried to use an expression build by me, but in the Definition Query, i am limited to a certain number of characters. I want to subset points from over 200 000 unique codes.

The codes are numeric. Here is a sample code:

IdRPL in (1607130,1607131,1609481,1609482, ....., 1709809)


texts: Can the qgis print composer memory limitation be overridden?

texts:
I'm trying to export large images from the QGIS 2.16.3 print composer. As far as I can tell, whenever the calculation of dpi and export size exceeds ~1615MB, QGIS will show the message:

  "Trying to create image #1( 21239x25360 @ 317dpi ) may result in a
  memory overflow. Please try a lower resolution or a smaller paper size"

And then not attempt the export. The same limitation appears to be applied when trying to do the export from the command line as well.

Is this a standard limitation, or is it based on system specs? (I have 32GB of system RAM). Is there is a way to override it and make QGIS attempt to output the image?


texts: I have a "park area" shapefile and a "lake" shapefile both consisting of polygons. There are some areas where the park area polygons do not entirely match the borders of the lake polygons due to data collection inconsistencies, even though the park area should match the lake outline completely. I am wondering if there is a way to properly extend the park area polygons in certain locations so they follow the outline of the lakes. I have attempted to use the Draw functions in ArcGIS for Desktop to no avail. I am not overly concerned with the accuracy of the park area polygons as this editing is for cartographic purposes, so if the park area ends up "inside" the water bodies it isn't an issue, but if there is a way to make the park area polygon automatically end at a lake border that would be preferable. Below is an image showing the issue I am having - The yellow background "sticks out" underneath the park area polygon at points it should not.


texts:
I am trying to merge 8 .asc files using merge option in QGIS 2.4 but when I do that it makes like contour of the file when I want to merge the files together.

.ASC Files:

Result after merging:


texts: PostGIS polygon clip inaccurate using ST_Clip and ST_Union from tiled raster


texts: Is it possible to convert a .shp file to .gml?

texts:
Is it possible to convert a shapefile to a gml file?
If so, how and using what tools?


texts:
I have downloaded a shape file of Chicago building footprints from here. Unzipped it and ran

ogr2ogr -f csv -dialect sqlite -sql "select AsGeoJSON(geometry) AS geom, * from 'geo_export_98941210-c995-4502-8ebb-93a09f29d615'" chicago_building_footprints.csv geo_export_98941210-c995-4502-8ebb-93a09f29d615.shp

Why do I get the following error? Have no issues with several other shape files downloaded from the same site. Using GDAL 2.2.3, released 2017/11/20 on Ubuntu 18.04.2.

  ERROR 1: Cannot create virtual table for layer
  'geo_export_98941210-c995-4502-8ebb-93a09f29d615' : CREATE VIRTUAL:
  invalid SQL statement : CREATE TABLE
  "geo_export_98941210-c995-4502-8ebb-93a09f29d615"("date_bldg_"
  DATE,"time_bldg_" VARCHAR(254),"bldg_condi" VARCHAR(254),"date_bldg_"
  DATE,"time_bldg_" VARCHAR(254),"date_bldg_" DATE,"time_bldg_"
  VARCHAR(254),"bldg_id" FLOAT,"bldg_name1" VARCHAR(254),"bldg_name2"
  VARCHAR(254),"bldg_sq_fo" FLOAT,"bldg_statu" VARCHAR(254),"cdb_city_i"
  VARCHAR(254),"comments" VARCHAR(254),"date_condi" DATE,"time_condi"
  VARCHAR(254),"create_use" VARCHAR(254),"date_demol" DATE,"time_demol"
  VARCHAR(254),"date_edit_" DATE,"time_edit_" VARCHAR(254),"edit_sourc"
  VARCHAR(254),"edit_useri" VARCHAR(254),"f_add1" FLOAT,"footprint_"
  VARCHAR(254),"harris_str" VARCHAR(254),"label_hous"
  VARCHAR(254),"no_of_unit" FLOAT,"no_stories" FLOAT,"non_standa"
  VARCHAR(254),"orig_bldg_" FLOAT,"pre_dir1" VARCHAR(254),"date_qc_da"
  DATE,"time_qc_da" VARCHAR(254),"qc_source" VARCHAR(254),"qc_userid"
  VARCHAR(254),"shape_area" FLOAT,"shape_len" FLOAT,"st_name1"
  VARCHAR(254),"st_type1" VARCHAR(254),"stories" FLOAT,"suf_dir1"
  VARCHAR(254),"t_add1" FLOAT,"unit_name" VARCHAR(254),"vacancy_st"
  VARCHAR(254),"x_coord" FLOAT,"y_coord" FLOAT,"year_built"
  FLOAT,"z_coord" FLOAT,OGR_STYLE VARCHAR HIDDEN,"GEOMETRY"
  BLOB_POLYGON_XY_4326, OGR_NATIVE_DATA VARCHAR, OGR_NATIVE_MEDIA_TYPE
  VARCHAR) ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(select
  AsGeoJSON(geometry) AS geom, * from
  'geo_export_98941210-c995-4502-8ebb-93a09f29d615'):   no such table:
  geo_export_98941210-c995-4502-8ebb-93a09f29d615


texts: Automate the selection of duplicate geometries with QGIS3

texts:
I have a layer of buffers calculated from watercourses. Each watercourse/buffer has a Strahler Order (0 to 7).
I need to delete overlaps between buffers order 0 and buffers order 1 to start, but choosing to keep the highest order (1), that's why I can't simply apply "remove duplicate geometries", as I did first on seperate layers of buffers with the same order)
Here is a summuary of the several steps I apply, mainly using models :

create one layer for buffers of each order
apply v.clean tool, for area and the cleaning option "rmarea" that cuts polygons according to overlaps.
Before :

After :

Like I now have duplicated geometries at intersections, I apply remove duplicates

I did these steps for each order, and now I would like to remove duplicates of the lowest order. I tried the following steps :

merge buffers order 0 and 1
apply v.clean with the same parameters as below
select duplicates using the following expression, but it doesn't work :

intersects(  $geometry,  geometry( get_feature('Arthabaska_Buffer_RH_CE_D_20m_ordres01','ID_O_Str_1',0)))

In the previous selection, select buffers order 0 and remove them

What's wrong with this expression, and/or is there an easier way to do ?


texts:
I would like to write an array from a single band image to a 3 band vrt where each band has the ColorInterpretation() tag of red, green, blue respectively.

My code below;

from osgeo import gdal
import os

#IM = "/path/to/image.tif"
IM = im

### read image ###
ds = gdal.Open(IM)
X = ds.RasterXSize
Y = ds.RasterYSize
band = ds.GetRasterBand(1).ReadAsArray()

### write to 3 bands ###
driver = gdal.GetDriverByName("vrt")
outPath = os.path.join(os.path.split(IM)[0], "test_image.vrt")
outIM = driver.Create(outPath, X, Y, 0, gdal.GDT_Int16)
for i in range(1, 4):
    outIM.AddBand()
    outIM.GetRasterBand(i).SetRasterColorInterpretation(2 + i)
    outIM.GetRasterBand(i).WriteArray(band)
    print outIM.GetRasterBand(i).ReadAsArray()
    print outIM.GetRasterBand(i).GetRasterColorInterpretation()
outIM = None

The output is:

[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ...,
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
3
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ...,
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
4
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ...,
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
5

So setting the ColorInterpretation() (unlike when using GTiff, see here), but all values are zero. Why?!

#### UPDATE ####

I have been running this through ipython and just noticed this error on the command line:

Error 1: Writing through VRTSourcedRasterBand is not supported


texts: Using QGIS plugin that's designed for use in interface within PyQGIS?

texts:
I need to use a plugin that's designed for use within the QGIS interface in PyQGIS.

For example https://github.com/dgoedkoop/joinmultiplelines works great but I want to run it within a script, so that I can pass each group of lines to it to create the entire network.

So something like

layer = iface.activeLayer()
selection = layer.getFeatures() #all features
#selection = layer.selectedFeatures() #only selected features
dpr = layer.dataProvider()

completed_list=[]
for f in selection:
    print (f.id())
    text=f['RoadMntnc']
    if text not in completed_list:
        completed_list.append(text)
        line_selection=f.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"RoadMntnc" = text' ) )
        joined_line=MultipleLinesPlugin(line_selection)
        #merge layers based on https://gis.stackexchange.com/questions/169523/how-to-merge-layers-into-a-single-qgsvectorlayer/169543
        vl = QgsVectorLayer("LineString", "temporary_lines", "memory")
        pr = vl.dataProvider()
        feats1 = joined_line.getFeatures()
        for feature in feats1:
          pr.addFeatures([feature])
    else:
        pass

vl.updateExtents()
newfile=r'c:\scratch\roads_joined.shp'
QgsVectorFileWriter.writeAsVectorFormat(v1,newfile,"utf-8",v.sourceCrs(),"ESRI Shapefile")

How could I structure the joined_line=MultipleLinesPlugin(line_selection) to use the plugin?


texts: Create polygon of buildings from Quickbird images in QGIS

texts:
I have Quickbird/Ikonos/GeoEye image and want to extract all the buildings from the image as an independent layer in QGIS.
How can i achieve that?


texts: Flow mapping visualization using a road layer

texts:
I am trying to illustrate travel using a flow map on a road layer. My data is:

city of origin
destination
number of travelers

How do I use a polygon city layer (or point) and a road layer, with line thickness representing volume of travelers?

I am using QGIS.


texts: Data Driven map - Atlas, How to automate the filter feature?

texts:
I have a Postgres data source that needs to be filtered from hundreds of records to a few records, each time a make a map.  Is there a way to automate that, so that the atlas makes map pages for each 'unique' value in the filter field.
To be specific, I have list of states/provinces in the dataset where a particular language is spoken in the country. For each new language map, currently I am filtering the data source list using the filter function on the 'language' field.
Can Python commands or the Atlas feature be used to say generate such  filtered language maps in bulk ?


texts: How to remove tiny polygon outlines in geoserver sld?

texts:
The data I am using a polygon feature layer, consisting multiple polygons and I want to show the same without polygon outlines. In ArcGIS I simply use outline colour as no colour. This produce the below map.

However the same does not work in GeoServer SLD. I tried both, without any  tag definition, as well as the following:

<Stroke>
     <CssParameter name="stroke">#000000</CssParameter>
     <CssParameter name="stroke-width">0</CssParameter>
     <CssParameter name="stroke-opacity">0</CssParameter>
</Stroke>

But both of them produce the following:

How can I get rid of the tiny white line (or may be subtle gaps). I can not find anything like polygon-gamma, used in CartoCSS for the same, in SLD documentation.