texts:
I am trying to build a network analyst in flex, with the closest facility.

A problem has occurred during the process of solving closest route to find the incident.
In ArcGIS for Desktop, there is no problem about the distance between incident and facility, and the route can be solved easily.

However, when I publish in ArcGIS for Server and I set as the input data in flex, the distance is limit.  The closest route is only solved when the distance between incident and facility less than 1 meters.

Is there any parameter that I may have missed when I published the service, or is there an alternative procedure to manage the maximum distance of the closest route?


texts:
I have few rasters, they are on WGS84 4326, I want to save them as image and then to upload the image and that it will appear on the correct location.
for some reson, whenever I try to export the rasters as image, and the I try to upload it as raster it uploaded it but it doesn't display.

I tried also to merge all the rasters but all I got was black line  like this-

I'm afraid it may happen because some  rasters overlap a little on the edges , which is not problematic for my project but  maybe problematic for QGIS?


texts: Merge polygons based on Field A or Field B

texts:
I have a parcel layer I need to dissolve based on Owner OR Postal Address.
I know I can always dissolve based on Owner AND Postal Address, but that is not what I need.

So if:

owners are different but postal address is the same -> merge
postal address is different but owners are the same -> merge

I tried dissolving based on owner, maintaining statistics for postal address, and then dissolving based on postal address, but this is not entirely accurate because the dissolve chooses the first value of postal address it finds among equal owners.


texts: Source for Esri sample data from 911 Hot Spot Analysis

texts: I was trying to replicate the Hot Spot Analysis (Getis-Ord Gi*) from Esri, but cannot find the shapefile, like the "911Count.shp", "911HotSpots.shp", etc. Can somebody tell me where I can download these shapefiles?


texts: Choosing appropriate projected coordinate system for Africa

texts:
I am using the NNjoin tool to compute the distance between points and a multilinestring. Both layers are in CRS 4326 and cover most of Africa.

I've understood that I need to reproject my layers to a projected coordinate system to be able to get the distance between the points in metres rather than degrees.

How can I choose which projected coordinate system to use?


texts: Problem Inserting feature with Geotools WFS

texts:
With this tutorial

http://docs.codehaus.org/display/GEOTOOLS/Data+Writing

When I try to insert a feature

try {
DataStore geoToolsClient = getGeoToolsClient();

Transaction transaction = new DefaultTransaction();
try {

FeatureWriter<SimpleFeatureType, SimpleFeature> featureWriter = geoToolsClient.getFeatureWriterAppend(typeName, transaction);

            SimpleFeature simpleFeature = featureWriter.next(); // The problem is here!
              simpleFeature.setAttributes(featureToInsert.getAttributes());

              featureWriter.write();

              transaction.commit();

            } catch (Exception ex) {
                transaction.rollback();
            throw ex;
            }

This bug occurs at featureWriter.next(); statement which is supposed to provide an empty SimpleFeature. Then, we have to fill and save it.

This attempt to create an empty SimpleFeature by giving a value to NULL nillable attributes and default non nillable attributes.

We have a column of SDO_GEOMETRY type is not nillable in the bind table.

When the "next ()" attempts to create a default value for this attribute, the following method is called to create a default instance:

org.geotools.data.DataUtilities # defaultValue (Class Type)

But then the attribute is "com.vividsolutions.jts.geom.Geometry" type and defaultValue method does not know how to create a default instance of this type. This exception is thrown :

throw new IllegalArgumentException(type + " is not supported by this method");

And I’m not able to insert a spatial feature in the database

Thanks for any ideas!

Jeremie


texts: Registering database to ArcGIS for Server?

texts:
I want to register an SDE geodatabase to ArcGIS for server 10.1 . I use SQl server 2008 R2 ,
ArcGIS 10.1 . I backup my database from server "a" and restor to server "B" then I used registering database in ArcGIS publisher .I successfully registered database (fig1), but when
i want to publish my map ,the warning raised " datasource is not registerd with the server and data will be copied to the server"(fig2) . How do i get rid of the warning ?


texts: Calculating the average of an attribute for a points within a pixel

texts:
I want to calculate the mean of an attribute of all points within each pixel.

For example, if there are two points within a pixel like in the lower left most corner shown, I want to calculate (attribute 1 of first point+attribute 1 of second point)/2

This is just one cluster of points, I have many more across the image and cannot do it manually.


texts: Symbolizing Raster using single Unique Value?

texts: I'm dealing with a raster in GRID format. Currently it is symbolized in Unique Values and contains 5 classes. I'm trying to consolidate these to one class. Is there a way to do this and keep a single Unique Value?


texts: Capitalizing The First Letter Of Each Word In Attribute Field QGIS

texts:
I have tried capitalizing the first letter of each string in my field/column called "COMMUNITY" by referencing this post:

Capitalizing The First Letter Of Each Word In Attribute Field

Particularly:     value = <COMMUNITY>.title()

I get an Expression is invalid output. What am I missing?


texts: Line at map scale that do not vanish when zoomed out? (QGIS 3.8)

texts:
I'm drawing rivers and roads as linestrings with width set in map units from an attribute.

When zooming out, of course they vanish quickly, because 10, 20 or 100m is less than a pixel at 1:100,000

How do I tell QGIS to always draw that line, maybe by changing from map scale to millimeters or by setting a scale-dependent minimum width?

For comparison, Google Maps always draws major roads (highways, etc.) even when the scale becomes such that they're the width of cities. I'm looking for that kind of effect.


texts: National grid and OpenLayers do not match

texts:
I have points in a txt file with lon, lat and elevation. The CRS is the National Grid of my country (metric system). However, when I add them in QGIS and then try to lay over them a map from OpenLayers plugin, the points are 50 meters off their true location. How to fix this?

The national grid is Balkan Zone 7.


texts: Select specific polygon intersecting another polygon

texts: I have lots of building polygons in my parcels. Each with individual attributes. I want to select for each parcel the building polygon with the highest value from one field. How can I accomplish this in ArcGIS Desktop.


texts: Add polygons to PostGres to be queried by Nominatim

texts:
I'm trying to build out a geocoder that also returns the correct legislative district information from this data:
http://www.census.gov/geo/maps-data/data/tiger.html

I've got OSM data imported into a PostGres database and I'm using Nominatim to query that information, but I cannot figure out how to lay the above information on top of that and return it via Nominatim. After reading a lot, I'm assuming it would make sense to layer the legislative district information as polygons.

I'm not looking for the polygon to be mapped as the usage will be strictly API. What type of geodata file do I need to download? And how would I go about getting it included in the queried elements of Nominatim?


texts: How to export a list of intersections of 2 datasets of boundries/polygons

texts:
I am using the following datasets:

NYC Zip Codes: https://data.cityofnewyork.us/Business/Zip-Code-Boundaries/i8iw-xf4u
NYC City Council Districts: https://data.cityofnewyork.us/City-Government/City-Council-Districts/yusd-j4xi/data

I want to generate 2 lists, the zip codes in or partially in each council district, and the council districts in or partially in each zip code.

So far I have tried this in Carto but I am open to other solutions.

I have found examples of intersecting a second layer of data with points, but nothing about intersecting 2 layers of data where both data sets are boundaries.

Is this possible?


texts:
How does buffering objects in a raster data set differ from buffering in a vector data set?
What are the steps on making a raster buffer?


texts: Getting meters instead of degrees from GRASS v.distance?

texts:
I'm trying to measure the distance between points to a coastline. I have a shapefile of the "High Seas" which I got from NOAA's website here.

http://www.nws.noaa.gov/geodata/catalog/wsom/html/marinezones.htm

I also have a list of points which represent buildings.

I've played with different solutions, but I think that GRASS's v.distance is the best solution. I've been able to get my data into GRASS and can plot the map.

The problem is that when I run the v.distance tool, I'm getting decimal degrees, but I really want either meters, miles, or kilometers.

So the v.distance function seems to be the best option.

However, I've been unable to get meters out of this data. In the help guide, for v.distance, I did find this sentence.

In lat-long locations v.distance gives distances (dist and to_along) in meters not in degrees calculated as geodesic distances on a sphere.

Its not exactly clear to me what this means. It sounds like I should set up the location with a EPSG code that uses lat long. I've used a handful to no avail. For example, I tried EPSG 3857, 3786, 4326, all to no avail.

As a side note, I did try QGIS's distance to hub that's part of the MMQGIS plugin. It appears that maps to the nearest centroid of a polygon. I'm trying to get the distance to the nearest pixel/border of a polygon. Or the nearest pixel/edge of a line.


texts: Adding Rules to Topology with Arcpy: getting ERROR 000714

texts:
I want to make a very basic python script, however I keep running into error 000714 (general/unknown script error).

The topology is created, the features are added to it, however the rule is not. I've tried different rules, same problem; I will want to add more rules and validation if this error is overcome.

I have tried doing the same process in modelbuilder but the I get error 999999, the syntax is the same when I export the code from it.

Anyone run into this problem or find solution? I am using ArcGIS 10.1.

Code (The script tool's only parameter is a Dataset):

import arcpy, sys
from arcpy import env

env.workspace = sys.argv[1]

topolo = "GBA_Topology"

transLn = "TRANSPORT_LINE"
gbaTransLn = "GBA_TRANSPORT_LINES"

if arcpy.Exists(topolo):
    arcpy.Delete_management(topolo)

arcpy.CreateTopology_management(env.workspace, topolo, "0.001")

arcpy.AddFeatureClassToTopology_management(topolo, transLn, "5", "1")
arcpy.AddFeatureClassToTopology_management(topolo, gbaTransLn, "5", "1")

arcpy.AddRuleToTopology_management(topolo, "Must Not Self-Intersect (Line)", transLn, "", "", "")


texts: Configurable Map Viewer: copy of Search widget fails to reference renamed id="mySelect2"

texts:
I made a copy of the Search widget to use on a different layer.  I renamed the copy Search.js file to Search2.js, Search.html to Search2.html and referenced it like so in config/viewer.js:

   search2: {
            include: true,
            id: 'search2',
            type: 'titlePane',
            canFloat: true,
            path: 'gis/dijit/Search2',
            title: 'Find Precinct',
            open: true,
            position: 3,
            options: {
                map: true,
                mapClickMode: true,
                mapRightClickMenu: true
            }
        },

I changed the id of the selection box in the js and html copies to "mySelectp."

I get an error in the browser:

  dojo/parser::parse() error Error: Tried to register widget with id==mySelect but that id is already registered

I can't figure out why the Search2 widget is seemingly still trying to reference the element with id==mySelect, instead of the renamed id==mySelectp.


texts: Dividing map into features with invisible divisor lines using QGIS?

texts: I'm working with QGIS 1.8.0. Is there a way to divide a map with a grid but make the divisor lines invisible? I have a shapefile as the one in the left and I want it to see as the one in the right but still having the divisions (with the 82 different features).


texts: Prevent GRASS from converting utf-8 to latin-1

texts:
I am opening two shapefiles in QGIS 2.18 and then perform a v.overlay from the GRASS toolbox.
The shapefile contains attributes with characters that need Utf-8 to display all the characters correctly. However, when I use v.overlay, the algorithm tries to convert all the characters to latin-1 which does not work. A lot of characters are converted to ?.

How do I tell GRASS, when used from QGIS, to use utf-8?


texts: Export raster in pieces to stitch together afterwards or other way to print a big map?

texts:
I used QGIS to put together an OSM XYZ Tile-layer and some WMS-layers. Now I want to make a print map out of it. When I try to export it from the print layout editor to PNG or PDF in the desired resolution and size, QGIS quickly eats up all memory of the system and then gets killed or crashes due to out of memory.

Before exporting, it warns me with a message that it will take something like 500 MB of memory. But when running it will very quickly eat up the 6 GB of (remaining) system memory, then the 8 GB of swap, and then gets terminated.

Is there an option to make QGIS work in pieces and stitch them together afterwards (probably with some other tool), or any other way to export a big raster image from the print composer?

It is currently QGIS 3.6.1 on Arch Linux.


texts: Calculate raster mode with r.mode - problem with decimals

texts:
I am working with a raster file that represents the percentage of vineyards by commune (e.g., municipality) in a region in France. The statistics were calculated for each commune based on surveys from the ministry of agriculture. My problem is that the statistics and the resulting raster were elaborated by a colleague and our boundaries for different communes do not match. In other words, some cells which have the label for one commune according to my commune polygon vector (and rasterized afterwards) have the label of other commune in according to his map. I though I could create a new raster in which a unique value is given for each commune using the command r.mode.

The most frequent value within each polygon could be extended to all the commune area.

I am using GRASS GIS. My problem is that the raster file has float values and r.mode requires integer values.I need to transform my float values into integers, but before I need to multiply my values by 1000000 or similar since all my raster values are smaller than 0.

I would like to do it in r.mapcalc, using an expression like:

r.mapcalc value_int@covariates = int(value@covariates * 1000000)

but this expresion gives me an error.

I may not have a lot of imagination, so if any solution can be done in R, I would also be really thankful.


texts: Implementing Helmert transforms between reference frames / coordinate systems?

texts:
I'm trying to implement a 14-parameter Helmert transform to convert between terrestrial reference frames, and having difficulties with it.

I have successfully followed the example in ITRF to GDA94 coordinate transformations, Dawson & Woods, 2010 (Appendix A); d28rz98at9flks.cloudfront.net/71433/71433.pdf.

Using the coordinate in the sample:

[-4052052.36775,4212836.04113,-2545105.10906](ITRF2005@2010.4559)

I use the following to normalise transform parameters:

const T = { x: tx + ṫx*δt, y: ty + ṫy*δt, z: tz + ṫz*δt };
const R = { x: rx + ṙx*δt, y: ry + ṙy*δt, z: rz + ṙz*δt };
const S = 1 + s + ṡ*δt;

Reference epoch is 1994.0 giving δt=(2010.4559-1994)=16.4559; resulting transform parameters are

T: { x: -0.04270,    y: -0.01706,    z:  0.02881 }
R: { x: -1.17163e-7, y: -1.01576e-7, z: -1.03731e-7 }
S-1: 1.1474e-8

Which match those given in the sample (rotation is negated as “Australia assumes the rotations to be of the coordinate axes”).

Applying the transform as follows:

const x2 = T.x + x1*S   - y1*R.z + z1*R.y;
const y2 = T.y + x1*R.z + y1*S   - z1*R.x;
const z2 = T.z - x1*R.y + y1*R.x + z1*S;

I get

[-4052051.76142,4212836.19454,-2545106.01462](GDA94@2010.456)

which matches (to within a mm) the sample transformation.

However, applying the same calculation to the Onsala observatory example from Proj.4 (github.com/OSGeo/proj.4/blob/2aaf53/test/gie/more_builtins.gie#L357), I don't get the correct result.

From coordinate:

[3370658.37800,711877.31400,5349787.08600](ITRF2000@2017.0)

with a reference epoch of 1988.0 giving δt=29; resulting transform parameters are

T: { x: -0.00714,    y:  0.00007,    z: -0.00383 }
R: { x: -1.73563e-8, y: -2.28347e-8, z:  4.31484e-9 }
S-1: 2.2400e-9

I get

[3370658.25318,711877.42306,5349787.15877](ITRF93@2017.0)

as opposed to the correct

[3370658.18890,711877.42370,5349787.12430](ITRF93@2017.0)

which is wrong by [0.06428,-0.00064,0.03447].

The NOAA HTDP agrees with the Proj.4 result, to 3dp.

Could anyone help with what I'm doing wrong?


texts: GeoServer altering raster values of multi-band raster WCS request

texts:
I'm hosting some NAIP imagery (4-band) through Geoserver and while WMS requests show true pixel values properly, if a WCS request is issued, the raster GeoServer returns has pixel values of nodata where 0 would be.

For instance, where a pixel in the original raster would have RGBI values of 38, 0, 2, 56, the returned raster has values of 38, nodata, 2, 56.

From the logs (Verbose), I can't see any sort of operation being performed (e.g. warp) that would result in an effectively new raster, so I'm not sure what's going on here.  No input transparent color is declared in the layer properties.  The request looks something like:

.../geoserver/wcs?service=WCS&version=1.0.0&request=GetCoverage&coverage=TEST_RASTER_NAIP&bbox=491592.0,3474409.0,492178.0,3474761.0&crs=EPSG:26912&format=geotiff&width=586&height=352

The result is also the same if I build it using resx and resy instead of width and height.


texts: How to get Name of IRaster object using ArcObjects?

texts: I have an IRaster object as output from a function. I need to be able to get the name of that raster. Is there a way that I can do that? I know that I can't mask an IRaster as an IRasterDataset. If I could, I would be able to get the Name property from that (I think). Thanks!


texts:
I have two images from two different time periods that are not registered (don't align) as shown below. I understand that image registration is possible using OTB's FineRegistration feature. However I'm a bit confused as to what some of the required arguments mean. Specifically the exploration radius (x and y) and metric radius (x and y). I have referred the function documentation but it's not immediately clear what these params do.
Could someone explain what these parameters are and their selection criteria?


texts: Efficiently relate two large shapefiles for use in Python script

texts:
I have a large National Land Use Code polygon shapefile for the San Diego area as well as a large point shapefile representing the center of each pixel from the Land Use Code raster for the same area.

The point shapefile is about 213,000 records and the polygon shapefile is about 170,000 records.  I'd like to spatially relate the points to the polygons so that the output will have each polygon and the points that fall within each polygon and will then use it in a Python script making use of a GIS library that can read and analyze shapefiles.

Since the two shapefiles I am joining are so large and the join type is a 1-to-many the resulting shapefile comes out to around 41gb when stored within a geodatabase, and won't even finish writing out if I don't store it in a geodatabase and just a folder.

Is there any way to spatially relate these two files in a more efficient way? I have attached a snapshot of the two files overlaid on top of each other.  Some polygons only have one relationship, while others have tens and hundreds of relationships.


texts: QGIS Server - WMS Request GetPrint fails with ATLAS_PK

texts:
I'm a junior developer and I'm working on a project trying to refresh an old GISWEB, to do so i chose QWC2+QGISServer.

Actual setup is:
Debian 10, Apache2, QGIS-Server 3.8, using .shp as layers... the QWC2 works fine, WMS works fine, i can print some templates, print png, etc...

The issue is that I'm trying to use the &ATLAS_PK (implemented in QGIS Server 3.6) and I'm getting an annoying "Internal Server Error". When I use "&ATLAS_PK=1" or "&ATLAS_PK=1, 2, 3", log says:

"CRITICAL Server: An error occurred during the Atlas print
WARNING: finish() called twice"

As I read in the code here it could be caused by 2 things :

"if ( pkIndexes.size() < 1 )
    {
      throw QgsException( QStringLiteral( "An error occurred during the Atlas print" ) );
    }"

or:

"QString errorString;
    atlas->setFilterExpression( filterString, errorString );
    if ( !errorString.isEmpty() )
    {
      throw QgsException( QStringLiteral( "An error occurred during the Atlas print" ) );
    }"

I can't determine which one is causing it.

If I try to use "&ATLAS_PK=*" it seems to work but I get "timed out" because I have like 4k features to print (I already changed the feature limit in QGIS Server conf). I'm going to try with a layer with less features as soon as I can and post updates.

GetProjectSettings WMS request returns that I have 2 templates both with atlas enabled.

There's not much documentation to read about, as it's a new feature...

UPDATE:
I tried to use atlas_pk with postgis data instead of .shp data and it works perfectly. if I have time I'll try to check if it's a problem with the .shp data or with the implementation of the atlas_pk that doesn't allow the use of .shp data.


texts: Geometry of edges created after using osm2pgsql

texts: I used osm2pgsql to convert osm data to postgresql. Now I have a table with edges and they have geometry like multilinestring. I am not sure what is meant be geometry of an edge. An edge is supposed to be a line between two points/nodes(start and end). How does the geometry thing come in? Can anyone please elaborate? And I want to know if routing functions like pgRouting actually use this geometry when calculating/displaying the shortest route?


texts: How to link R script to Algotithm output in QGIS modeler?

texts:
I constructed a simple R script for QGIS modeler (.rsx) that calculates groupwise summary statistics of a vector layer table (sum by field "subject" and "FRC").
The script correctly works. When executed, a windows is opened and you can manually choose the vector layer and the other parameters.
However, if I insert this script in the graphical modeler I am able to choose the output of a previous Algorithm as input to my script, nevertheless I cannot select fields from its attribute table, as requested by my script.

How can I do that?

----Edit-----

As an example, consider I want to construct a model that calculates the % of different land use types inside polygon buffers constructed around some street addresses.

The model represented in the picture below intersect the buffers with a polygon layer representing land use areas, then calculates the area of each new feature in the intersection layer.

Now I want to link my R-script to the output layer "Intersection_with_area" to summarize the total area for each address and for each land use type. Here below the .rsx code:

##intersection=vector
##address=field intersection
##soiltype=field intersection
##area=field intersection
##folder=string
##filename=string

library(doBy)
setwd(as.character(folder))

inters_df=as.data.frame(intersection)

summ=summaryBy(list(area,c(address,soiltype)),FUN=sum, data=inters_df, keep.names=T)
write.csv(summ,as.character(nomefile))

I am unable to select the fields "address", "area" and "soiltype" from the vector layer "Intersection_with_area" that is in my model.


texts: How to show attribute table of object when user clicks on it?

texts:
sorry if I'm asking stupid questions but I really spent days trying to figure out how to add functionality to my map so that when user clicks on some map object he gets the attribute table on the side with the attribute data of that object. Is it even possible? I don't speak about popups, I speak about actual table.
This is how my map looks now:

I would like that space on the east side of monitor populate with that attribute table when users click. I'm not saying that's the best solution but that's how I imagined that. But I'm open for suggestions.

As you can see I use geoext, openlayers and extjs and honestly, I don't know lot about php or some other frameworks so I would be happy if the solution could be solved using geoext or smilair. I was looking around but nothing what I found is what I'm actually looking for.

One more time, sorry if my question is unprofessional or something but I'm really kind of desperate :)


texts: How to install mapserver on Windows/Apache?

texts:
We need to get an internal instance of MapServer up and running. We want to use Apache and not IIS because we already have an Apache 2.4 instance up and running on our office intranet.

The Windows builds are available from a number of websites on this page: http://mapserver.org/download.html#windows

I chose the Gisinternals build because it has daily builds, and the packages are available as Zip distributions.

I need to run MapServer as a CGI from Apache.

There is a lack of accurate information on the internet about how to get this installed, configured and running correctly with Apache on Windows and I have been struggling for the last 2 days without success. Searched for over 6 hours today, trying various configurations and still got nowhere. I'm very exhausted frustrated.

If you have managed to get MapServer up and running with Apache on Windows, please tell me how you did it, which build you used and step by step configuration details. I also want to use it with mod-fastcgid.

This could become a great community post that would help others in the same situation.

UPDATE:
I got it up and running by copying the DLL files to the same CGI-Bin directory where the mapserver executable resides. However now I get the error:

  mapserv(): Web application error. Traditional BROWSE mode requires a
  TEMPLATE in the WEB section, but none was provided.

The url I browsed to looks like this:
http://localhost/cgi-bin/mapserv.fcgi?map=c:\mapservproj\proj1.map

Still looking for a good step-by-step guide to setting this up though..


texts: How can I carry over attributes from a GIS .shp to a SnakeGrid MicroStation .dgn, using FME?

texts:
I have the SnakeGrid projector so can convert coordinate systems but can't seem to carry the attributes over (I have tried AttribueManager and AttributeCopier).

Also, can I input several .shp and run a batch to export as .dgn with attributes carried over?


texts: Plot a range of values of a NDVI calculated in Google Earth Engine

texts:
I just perform a NDVI classification, and I want to know if I can plot or display the values of the NDVI lower than 0.8, or within a range for example from 0.3 to 0.5 and the rest of values in black or just not to plot them.

Here is part of my code:

//cargar coleccion de landsat
///////////////////////////////////////////////////////////////
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')

//geometry with the AOI
var roi = geometry;

//filtrar la imagen de acuerdo a la cobertura de nubes en el poligono
////////////////////////////////////////////////////////////////////
var l8f = l8.filterBounds(roi);

var withCloudiness = l8f.map(function(img_zone) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(img_zone).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: roi,
    scale: 30,
  });
  return img_zone.set(cloudiness);
});

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));

var image = ee.Image(filteredCollection.filterBounds(roi)
    .filterDate('2017-03-30', '2017-07-30')
    .sort('CLOUD_COVER')
    .first());

var red = image.select('B4');

//infrarrojo cercano
var NIR = image.select('B5')

var nvdi = (NIR.subtract(red)).divide(NIR.add(red)).rename('NDVI');

Map.addLayer(nvdi,{min: -1, max: 1},'ndvi_rice_1')


texts: Find the theoretically fastest route based on many partially overlapping GPX tracks

texts:
Between the points A and Z in a city, there are several possible routes that may partially overlap or cross at crossroads (such as A→B→C→D→Z, A→J→C→D→Z, A→K→C→L→Z, A→M→N→L→Z, etc.). This may be represented by a directed graph.

I have a lot of GPX tracks (timestamp, lon, lat with GPS accuracy) from A to Z following any of these routes. Some of them follow the same path, but may have different durations at different sections (because of traffic, lights, weather, works, etc.).

I'd like to find the fastest (less time) route between A and Z. It is possible, that this route has never been traveled by any GPX track entirely, but consists of edges that have been traveled by different routes.

What tools would you use to dissect any number of GPX tracks into sections given a list of nodes? (open source and/or some Python modules preferred)

Let's assume these nodes are at least 300 meters from each other and the dissection should happen exactly in the middle between entering for the first time and leaving for the last time of the buffer of 50 meters around the given point. Or any other deterministic algorithm.

The result could be in the form {('A', 'B'): [50.1, 55.0, 61.1, …], ('A', 'J'): [81.1, 101.7], …}, with a list of durations in seconds between two nodes. Then I'd use Python's networkx to find the theoretically fastest route.


texts: Building QGIS with Visual Studio 2015 on windows x64?

texts:
I tried to follow the guide at github, hence I encounter problem after problem (Now I have problems with: MSBuild.exe  to get the value of VCTargetsPath:, earlier I couldn't find any c++ paths etc.) I assume that I have done some mistake and missed something running the installation.

I'm curious if there is anyone that recently have build QGIS on Windows(10) with Visual Studio(2015) and can give some advise what to change in the .bat file, where I can set the path to bison and flex in cmake, what and where I shall use the Ninja.exe file that I downloaded according the guide.


texts: Data entry and manipulation using web pages

texts:
I am working on a project where we want to move a database from Microsoft Access to PostGIS, and we want the data entry forms to be web page based rather than forms within the Access database.

I am in the process of migrating the database from Access to PostGIS and I have started to use Geoserver, OpenLayers and GeoExt to create maps on pages to look at the data. All this is working fine.

My question is what are the best approaches to allow users to query/edit/enter the attribute and spatial data on the web page? Does anyone have any good links to any tutorials?

Any help would be much appreciated.

Many thanks,

James


texts: How to force extent of multiband image to match the extent of another image?

texts: I have two 8bit unsigned integer 4-band raster images that overlap (Raster 1 and Raster 2). One of the images is slightly smaller than the other (Raster 1). How can I extend the extent of the smaller raster (Raster 1) to match the extent of the larger raster (Raster 2)? The extended area in Raster 1 should contain 0 values for all the bands. NoData values are not wanted in this case. I use ArcGIS/Spatial Analyst, QGIS and any open source tools.


texts: does someone have experience in fixing geometry errors in QGIS 3.10 LTR? I'd like to fix geometries first and then use processing algorithm "Snap geometries to layer". So far, the algorithm doesn't work as there are still geometry errors. For fixing the geometry errors I used the processing algorithm "Fix geometries".


texts: arcpy PolygonToRaster ERROR 000918: Cannot retrieve feature class extent

texts:
I am trying to convert a polygon feature class (created intermediately as part of a larger process) to a raster, which will have the same snap settings and cell size as another input raster.  However, every time i try to run it, I get this error:

ExecuteError: ERROR 000918: Cannot retrieve feature class extent.

I cannot find any help on this issue and I don't know what to do.  I have already seen this question and it did not help resolve the issue.

This is the code:

cell_size = float(GetRasterProperties(in_raster, 'CELLSIZEX').getOutput(0))

env.snapRaster = in_raster

PolygonToRaster(buffer_fc, z_field_name, stamp_raster, cellsize=cell_size)

return None

The same parameters work when I execute "Polygon To Raster" diretly from the ArcToolbox GUI.


texts: I need to calculate regional slope for a coastline using a ~500 meter DEM covering a coastal plain and continental shelf where the slope for each grid cell is calculated by defining elevation extremes within a 10 km radius for each individual grid cell. Normally slope is calculated using a 3x3 window. How do I increase that window? Is there a way to do this in ArcMap?


texts: I have a point layer in PostGIS containing point locations of vehicles acquired by GPS over a prolonged period of time. The layer contains about 600k points. When loaded into QGIS points of this layer that are expected to be grouped along roads centerline are visualized as regularly spaced in a kind of grid. It looks like this is only QGIS visualization, these are not true points locations. This might be a well-known feature, but I failed to google anything similar. Can this be turned off and true point locations shown in QGIS?


texts:
We are plannig to use QGIS in our company to reduce software coasts.

Do we have to take permission or make a payment?


texts: Saving levelplot to a TIFF raster without changing pixel size?

texts:
I am working on change detection and have been able to visualize the results that I have created through levelplot. Previously I used to use Writeraster for storing my plots as tif files which did not manipulated my data in terms of pixel width or height. With levelplot, I am unable to use this function,as it throws error missing trellis, character function.

How can I save a levelpot as a TIFF imagery without legends to be opened in QGIS so that I can overlay it over basemaps?

Here is how I am creating levelplot:

changeDet1 <- calc(stack(lc1,lc2), fun = change)
codes_ <- data.frame(ID = grid_$code,value = paste0('from ',grid_[,1],' to ',grid_[,2]))
logical_test <- which(grid_$change == T) # remove no change classes
codes_ <- codes_[logical_test,]
codes_
a<-subset(codes_, ID %in% c(2,10,12,13,14,15,16,17,18,20,29,38,47,56,65,74))
rat <- levels(changeDet1)[[1]]
rat[["Changes"]] <- codes_
levels(changeDet1) <- rat
levelplot(changeDet1, par.settings=PuOrTheme(), xlab="", ylab="")


texts:
I am implementing a C++ application that must produce DEMs from LIDAR data. I am following the next workflow:

Write .las files from lidar data.
Use PCL algorithms to transform/filter the point cloud.
Use GDAL/OGR to create the raster/grid.
Create a 3D model (.obj) from raster.

In this process I am doing the next file conversions:

From ".las" to ".pcd" (To work with PCL).
From ".pcd" to "vector format" (i.e. "ESRI Shapefile") (To work with GDAL).
From "vector format" to "raster format" (i.e. ".tif") (To generate the DEM).
From "raster format" to ".obj" (To visualize the DEM in 3D).

I'm using gdal_grid to create the grid from vector format file and gdal_dem to create the DEM using GDAL interpolation algorithms.

I'm not sure If this is the best way to create a DEM from .las files. I want to use third-party libraries to transform, filter or interpolate the point cloud and PCL/GDAL have these tools.

Anyone with experience with DEMs can give me information / advice?


texts: Offsetting symbology position in ArcGIS Desktop?

texts:
My data of polylines has lots of attributes to be rendered with a symbology, I made it such that each attribute of interest would be a layer. I would like to display each of these layers at the same time so that the attributes could be seen at the same time. Since the top layers cover the lower ones, I set the symbology so that the layers would have increasing different thickness, with the most bottom layer having the thickest line.

Is there a better way to do this, like offsetting a layer's symbology a little bit to the left or right so that each layer can be displayed side by side? Visually, something like this.

I'm using ArcGIS 10.


texts: Updating the value of a field with pre-defined values

texts:
I am trying change the value of a row for a field whose type is "short" in arcmap using a python script.

with arcpy.da.UpdateCursor(FILENAME, FIELD) as cursor:

    for row in cursor:
        row = 'Center'
        cursor.updateRow(row)

with this i receive the error "RuntimeError: The value type is incompatible with the field type."

I know I am giving a string type to a "short" type field. Is there a way around this? There are pre-set values in this field. Do I need to run some other kind of code to create a list of possible values for this field, then select the value based on some kind of list interaction eg row = fieldlist[1]

Please ask me to clarify anything, but this is everything that seems relevant to me at this point.


texts: Tassled Cap Transformation of Landsat-8 in ArcGIS using surface reflectance output from the USGS

texts:
I see lots of references for Tasseled Cap transformation for Landsat-8 Top of Atmosphere reflectance.  What I just got from Earth Explorer is a product called "surface reflectance", which sounds like it's atmospherically corrected and not the same product as the "top of atmosphere reflectance" product.

Should I re-order a different product?

Or, are there updated coefficients I might use for what I have?

Or (shudder), is there a tutorial on how the develop my own coefficients?


texts:
I got a DXF file with some height data which I can easly open in QGIS, the problem is when I import it to QGIS is that the "float" is separated from it's Integer. Instead of one point with 4.8 is there two points one with 4 and one with 8. And those points that is 5.0 does not have its "float" and is just "5". Like this:

Is there any tool in QGIS that I can use to combine the float and integer? The "float" is approximativ 1.7m south east of the integer.


texts:
I'm using overpass to retrieve information about OSM ways around some coordinates.

[out:json][timeout:25];
way["highway"](around:100,49.8385,24.0199);
out meta;

The default meta output is not sufficient for me -- I need to have city, country etc for every returned object.

I have two ideas on how to solve it:

Use recursion up (<) to traverse all the objects up
Use is_in query inside a foreach command for each object

I'm having problems with syntax / implementation on Overpass QL for both the ways.

My code for #1: (https://overpass-turbo.eu/s/NUi)

[out:json][timeout:25];
(way["highway"](around:100,49.8385,24.0199);
<;
);
out meta;

The result looks like recursion isn't applied.

My code for #2: (https://overpass-turbo.eu/s/NUj)

[out:json][timeout:25];
way["highway"](around:10,49.8385,24.0199);
foreach(is_in; out;)

Could you please point me on what is wrong here?


texts:
I'm looking for the best coordinate system projection to use for cartesian calculations while a vehicle is running from its starting location to about 100 miles from that spot.

The offline format for coordinates will be WGS84 geodetic coordinates. For online operation I am looking for the best coordinate system projection to use. We have decided that the best way to go about this is to define a projection with the origin at the vehicle's starting location.

Currently the primary candidates are an orthographic aka standard cartesian coordinate system set tangent to the earth at the starting lat/lon of the robot, or a transverse mercator projection defined at the starting x/y coordinate of the robot. We are avoiding UTM coordinates because we are likely to be on or near UTM boundaries so distortions would be a problem.

I am leaning towards the transverse mercator coordinate system, though I am not 100% confident how to appropriately configure the datums, particularly the scale factor, at arbitrary lat/lon coordinates.

For orthographic it seems that distortion would become a problem as early as the 50 mile mark from the starting location.

Ultimately my question is which coordinate system would be the best choice, and how to I choose parameters such as the scale factor appropriately based on the starting location?


texts: GeoServer - change Master Password - masterpw.info missing

texts:
I'm using GeoServer 2.5.1 and I want to change my default Master Password as is suggested by GeoServer itself when I log in:

  The master password for this server has not been changed from the default. It is highly recommended that you change it now. Change it

The Change it link leads me to the Change password form, though I need the default Master Password for this. According to this doc, I'm supposed to find it it at security/masterpw.info, however the only files I find there are masterpw.digest and masterpw.xml.
I've tried a bunch of sequences that looked like pws from these files and a file passwd in security/masterpw/default but all I am ever getting is

  org.geoserver.security.validation.MasterPasswordChangeException: Current master password invalid

This is reproducible on two different servers and plenty of fresh installations.
Am I overseeing something? How can I change my masterpw without masterpw.info?


texts:
After some package updates on my Debian machine all my database queries involving postgis functions or types are returning:

Error message from server: ERROR:  could not load library "/usr/lib/postgresql/9.1/lib/postgis-1.5.so": /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.18' not found (required by /usr/lib/x86_64-linux-gnu/libstdc++.so.6)

Is there a way to resolve this issue? I have to move spatial data between machines, but pg_dump also fails. Is there another way to safely dump my postgis databases?


texts: Fixing tiny microscopic Map font size on UHD?

texts:
QGIS 3.0
Added OSM layer via QuickMap services. This is what I see on my UHD display:

I understand that tile server just seeds tiles and we don't have any control over tiles generation. But why not to scale things if you're on UHD display?

I would expect QGIS to handle this by itself, we're at the 2nd quarter of 2019 after all.

Any ideas how to get it working?

It's impossible to work with maps in such circumstances.


texts: IfeatureLayer.search() second time return error

texts:
I have the following method in my code and have a button that when clicked run this code. In the first click all of code runs without error. but in second run I receive an error on line:
IFeatureCursor fcDistrict = flDistrict.Search(null, false);

Error:Exception from HRESULT: 0x8005018B

Void function(){
    IFeatureLayer flCircuit = ClsGeneral.FindLayer(ClsShares.FC_CIRCUIT);
    IFeatureLayer flSubstat = ClsGeneral.FindLayer(ClsShares.FC_SUBSTAT);
    IFeatureLayer flLoad_center = ClsGeneral.FindLayer(ClsShares.FC_LOAD_CENTER);
    IFeatureLayer flDistrict = ClsGeneral.FindLayer(ClsShares.FC_DISTRICT);

    IFeatureCursor fcDistrict = flDistrict.Search(null, false);

    IFeature fDistrict = fcDistrict.NextFeature();
    ISpatialFilter sf = new SpatialFilterClass();
    while (fDistrict != null)
    {
        int indCFCapacity = fDistrict.Fields.FindField("HIGHDIST_CIRCUIT_FREE_CAPACITY");
        int indHSFCapacity = fDistrict.Fields.FindField("HIGHDIST_SUBSTAT_FREE_CAPACITY");
        int indTSFCapacity = fDistrict.Fields.FindField("TRANSFER_SUBSTAT_FREE_CAPACITY");
        fDistrict.set_Value(indCFCapacity, 0);
        fDistrict.set_Value(indHSFCapacity, 0);
        fDistrict.set_Value(indTSFCapacity, 0);

        sf.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
        sf.Geometry = fDistrict.Shape;
        sf.GeometryField = flDistrict.FeatureClass.ShapeFieldName;
        IFeatureCursor fcCircuit = flCircuit.Search(sf, true);
        IFeature fCircuit = fcCircuit.NextFeature();
        while (fCircuit != null)
        {
            int indVol = fCircuit.Fields.FindField("UTNOM_VOLT");
            int indFCCapacity = fCircuit.Fields.FindField("FREE_CAPACITY");
            int Vol;
            int.TryParse(fCircuit.get_Value(indVol).ToString(), out Vol);
            if (Vol < 132)
            {
                double oldValue;
                double.TryParse(fDistrict.get_Value(indCFCapacity).ToString(), out oldValue);
                double sumValue;
                double.TryParse(fCircuit.get_Value(indFCCapacity).ToString(), out sumValue);
                fDistrict.set_Value(indCFCapacity, oldValue + sumValue);
            }
            fCircuit = fcCircuit.NextFeature();
            }
        }
    }
}


texts: Trying to map through and reduce regions of image using Google Earth Engine?

texts:
I'm working on building a model in google earth engine that will estimate carrying capacities for bird species based off of a table I have of bird densities per land cover type and an image of land cover types.

I'm at a point in my process where I have an image of bird densities across my region of interest. I want to use reduceRegion() with an ee.Reducer.sum() to get the total carrying capacity across my region of interest. My code looks like this:

//Density_Image is the image of bird densities I want to reduce.
//ROI is a FeatureCollection of regions within my region of interest.
print(Density_Image.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: ROI.geometry().geometries().get(3),
    maxPixels: 1e11,
    tileScale: 16}));

This code works when I feed in a single feature (~1e8 pixels), but fails when I try to feed in the area covered by my entire feature collection (~2e9 pixels) because the computation takes too long.

To get around this issue, and also because I would like to store the carrying capacity of each subregion, I'm trying to write a .map() function that will loop through each feature (region) in my ROI FeatureCollection, calculate carrying capacity for that feature, and then store carrying capacity in a new property I'll set within that feature, 'pop_est'. Here is my code for that:

var pop_est_table = table.map(function(f) {
    print(f);
    // line above throws an error that argref = null.
    Density_Image = ee.Image(remap2);
    print(Density_Image);
    var pop_est = Density_Image.reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: f.geometry(),
        bestTry: true,
        maxPixels: 1e12,
        tileScale: 16});
    return ee.Feature(f).set('pop_est', pop_est);
});

The code I've written has at least two things wrong with it.

1) The code will run and set a new property in each feature class called 'pop_est', but these values are always 0. Additionally, trying to print f throws an error:

  Feature (Error)
  Failed to decode JSON.
  Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
  Object: {"type":"ArgumentRef","value":null}.

I believe this error occurs because I am trying to call print(), a client-side function inside a server-side function, map(). I have a hunch the first error is also a result of some client-side/server-side confusion and would love some help on how to fix my code.

2) Even if this code does work I am worried that calling a reduceRegion() on each feature will take too long, causing google earth engine to shut down the computation.

How do I allow for a longer runtime or devote more computing power to this operation?

I've read that exporting is one way to attain longer runtimes, but I'm unfamiliar with the procedure.

I've tried several other approaches at this point, like clipping each feature out of my Density_Image and then putting those in an ImageCollection and mapping through that. I also tried reduceRegions(), but neither worked because their computations took too long.


texts:
I have a set of vector lines in grass which I need to convert into polygons.  However, because of bad base data, the polygons don't always close.  To fix this, I need to extend each disjointed line by a certain distance and connect it to the first line it hits.  The problem is that I can't figure out how to do this.  I have experimented with v.clean, v.generalize, v.edit, and (converting to a thinned raster) r.grow, r.neighbors and even r.flow, but I can't find anything which will let me do this.

Can anyone point me to a tool that would do this?


texts: When does ArcGIS Desktop set the metadata sync attribute?

texts:
I have some metadata that is being managed by the standard ArcGIS Metadata tooling accessed via the description tab in ArcCatalog.

Unfortunately (at least) one of the fields is not being applied correctly and when I look at the underlying xml I find that the element has its sync attribute set to false:

<formatName Sync="FALSE"/>

I remember that when we were managing our metadata via custom tooling we would have to be mindful of the esriMetadataSyncAction passed to the IMetadata.Synchronize() method or we would find the updates were only partially applied.

This leads me to ask, under what circumstances does ArcGIS set the sync attribute? And is there any way to override this using the standard tooling?


texts:
Is there a way batch remove or delete null geometries from hundreds of small polygon vector files in QGIS, while saving the new clean file to their original directories? I can't check for and fix geometry errors individually but I'm happy to give up some precision / data integrity / knowledge for a method that saves time.

These shapefiles will ultimately need to be processed but any bulk batch processing I do throws up null geometry errors. Until I fix these there's very little I can do.

I'm wondering whether 'saving as' a different file type might remove the error? If so, is there a way of pointing each of the output files back to their original directories without having to do this for each one?


texts: Map output with multiple "Zoom Levels" in QGIS 3

texts:
Just wondered if it is possible to do BBC-style multiple "Zoom Levels" output in QGIS 3 like this:

I used to use another app such as Adobe Illustrator to overlay the maps of different scales. However, as the maps were independent, it was neither fast (another step in another App) or precise (centres' coordinates may differ).

In the report I am working on, I intend to analyse a number of cases around the world, each of which is accompanied by a group of maps like the above example.

How can I have one underlying map for each location and export, say, three maps of different scales?


texts: Python standalone on linux: No module named qgis.core

texts:
I am using the Spyder IDE on my Ubuntu LTS 16.04. It does not find the module qgis.core. Therefore I wanted to manually add it to its search path but I do not now where it is located. What is the default path to that module on ubuntu?
I am using QGis 2.16

EDIT:
I tried to add this path:
/usr/share/qgis/python but it did not help

EDIT_2:
I've printed out my sys.path to show me all paths considered by the interpreter:

['', '/usr/share/qgis/python', '/usr/lib/python2.7/dist-packages/qgis/core', '/usr/share/qgis/python/plugins/processing/core', '/home/matt/anaconda2/lib/python27.zip', '/home/matt/anaconda2/lib/python2.7', '/home/matt/anaconda2/lib/python2.7/plat-linux2', '/home/matt/anaconda2/lib/python2.7/lib-tk', '/home/matt/anaconda2/lib/python2.7/lib-old', '/home/matt/anaconda2/lib/python2.7/lib-dynload', '/home/matt/anaconda2/lib/python2.7/site-packages/Sphinx-1.4.1-py2.7.egg', '/home/matt/anaconda2/lib/python2.7/site-packages/setuptools-23.0.0-py2.7.egg', '/home/matt/anaconda2/lib/python2.7/site-packages', '/home/matt/anaconda2/lib/python2.7/site-packages/IPython/extensions', '/home/matt/.ipython']

Do I miss something important there?


texts: Overlaid raster gets shifted around when printing to pdf or when resizing graphics device with mouse

texts:
I've spent the past couple of days learning how to use the maps, raster, and rgdal packages to create nice images, but I'm running into trouble when trying to make pdf's for actually presenting.  I'm overlaying two rasters -- one has my data and the other is a mask.  The problem is that the mask moves around depending on the size of the output pdf, AND also with simple resizings of the graphcis device.

I've made a little minimal example here:

library(raster)
library(maps)
library(mapdata)
library(rgdal)
lat = seq(from=36,to=40,by = .01)
lon = seq(from=-121,to=-117,by = .01)
d = expand.grid(lon,lat)
colnames(d) = c('lon','lat')
d$z = sin(d$lat)-cos(d$lon)+rnorm(nrow(d))
head(d)
coordinates(d)=~lon+lat
gridded(d)= TRUE
r = raster(d)
dd = expand.grid(lon[lon<120],lat[lat<37])
colnames(dd) = c('lon','lat')
dd$z = 1
coordinates(dd)=~lon+lat
gridded(dd)= TRUE
rr = raster(dd)
map('usa',xlim=c(-123,-116),ylim=c(35,41))
plot(r,add=TRUE)
plot(rr,add=TRUE,col="black",legend = F)
dev.copy2pdf(file = "example.pdf",height=8,width=8)

Here is the output when first rendered on my screen:

And here is the output when I resize the window!

and the example PDF is different again!  (Unfortunately SE doesn't seem to support pdf images, but if you run the code you can see it on your own machine)

Any idea how to work around this behavior?

Thanks!


texts:
I am new to python and I am writing my first plugin that add selected layers from a listwidget to QGIS.
So good so far, but it does not allow me to select several layers at once. Am I missing something in my code for that option?

 def run(self):
    # First clear the list widget
      self.dlg.listWidget.clear()

    # Fill the list
      f = codecs.open("C:/temp/testlist.txt", 'r', 'utf-8')
      for line in f: # Repeat for every line in the list file
        fileData = line.strip() # Remove any new line characters from the line (strip)
        if fileData != "": # Ignore empty lines
            if fileData[:1] != "#": # Ignore lines that start with hasch-tag
              tempList = fileData.split(",") # Separate the line by comma signs and put it in a temporary list
              layerProvider.append( tempList[0] ) # Add the first item to the layerProvider list
              layerName.append( tempList[1] ) # Add the second item to the layerName list
              layerSource.append( tempList[2] ) # Add the third item to the layerSource list
              self.dlg.listWidget.addItem(tempList[1]) # Add the Name of the layer to the list widget.

      f.close()

    # show the dialog
      self.dlg.show()
        # Run the dialog event loop
      result = self.dlg.exec_()
        # See if OK was pressed
      if result == 1:
            # do something useful (delete the line containing pass and
            # substitute with your code)
           selectedLayers = self.dlg.listWidget.selectedItems()


texts: Getting GPSTool icon back on manage layers tool bar of QGIS?

texts:
The GPS Tools icon has vanished from the "Manage Layers" tool bar and I would like to get it back.

I have tried removing and re adding the tool bar but that did not work and I know that I can access it from Vector -> GPS (but I keep forgetting ;)

Hmm... I think it must have gone when I upgraded to QGIS 2.18.2


texts: Calculating density map of cycleways using fishnet in ArcGIS Desktop?

texts:
How can I calculate density map of cycleways (line feature) which refers to net grid in ArcGIS?

Result by using function Line Density in ArcGIS for me is not satisyfing.

I created reference grid net by using Fish net. Now I would like to know how many m or km of cycleways occure in every grid to create density map.
I enclose youtube link from QGIS program with result which I would like to achieve in ArcGIS.

https://www.youtube.com/watch?v=-nb4-1ZRJbE

Is it possible?


texts: How to credit layers fom OpenLayers plugin for pdf/printed maps?

texts:
Anyone who has some insight into copyright law, went already through the terribly long terms of use and able to describe how should be the layers from OpenLayers plugin referenced / credited in PDF or printed maps? In the plugin there are links as follows:

The terms of service of the used map providers apply!

Selected links:
OpenStreetMap,
Google Maps,
Microsoft® Bing™ Maps

But after scan through these it all appears to me not saying anything about use for printed maps or PDFs. I must admit though I haven't read the whole thing word by word yet and hope I won't have to.


texts:
I have been having a conversation at the office about whether to use ground control when using drones/UAVs to capture aerial images. The ground control that I am referring too is any permanent survey markers, or Real Time kinematic systems that would aid in improving the spatial accuracy of the data captured. I would imagine that the level of ground control you would use, would depend on the accuracy that you were trying to achieve.

Would anyone know if there are some guidelines or provide some advice, for the use of ground control within Australia, when conducting an aerial surveys (UAVs) for use within photogrammetry?


texts: Can anybody tell me how I can combine the geo objects of one relation to one single object? I reckon the Object Combine function is right one, but I have no idea how to specify the code...


texts: Has ArcMap's background progress bar been fixed in 10.2?

texts:
I don't know if any other geoprocessing tool authors are pondering whether to invest the time and effort in a workaround for this "feature", but

in 10.0 background processing tasks had a progress bar indication both at the bottom of the main window and as a little 'clock' in the results window
in 10.1 the progress indication is gone, for background processes

Has anyone installed 10.2 yet?  Have we got it back?


texts: Is it possible to use PostgreSQL Geometric Operators for Latitude, Longitude?

texts:
I'm capturing tweets and saving them into a PostgreSQL table like:

CREATE TABLE tweets (
  id bigint,
  location point
)

And I'd like to select all tweets within a radius from a city, ie:

Buenos Aires (lat: -34.603722, lon: -58.381592, radius: 2Km)

As far I could understand I need to run a query similar to:

SELECT * FROM tweets WHERE location <@ circle '((-34.603722, -58.381592), 2000)'

I'm not sure if this approach makes sense, can I use Latitude/Longitude on simple PostgreSQL? Or I need to install PostGIS, I rather prefer not having more dependencies if it is possible.

I saw a few posts where people says a point should be (longitude, latitude) and I'm not sure unit is metric (meters) or imperial (miles).

Reference: Postgres Contained in operator


texts:
I have vector/topo/image tiled maps for some places in North America which were published as ArcGIS Server REST service. These maps (A) are based in EPSG:26912 projection. Now I want to use these as base map for my OpenLayers web app; besides, I also want to consume other base map sources (B), e.g., OpenStreetmap, Google Maps, ESRI, etc. Since most of these maps are in EPSG:900913 which is spherical Mercator projection to my knowledge, if I just add these base map layers into OpenLayers map component, the group A and group B won't overlay correctly, which means they are displayed as separate maps.

I know this is something related to projection, and might need some code with proj4js, however, as a newbie, I don't really know what to start with.

Basically, I want to put OSM/Google maps as the base map for the whole world in tier 1, then put EPSG:26912 base maps for some places in North America in tier 2, and put other WMS or WFS in the top tier.

do I need to transform EPSG:26912 to EPSG:900913 or opposite? I guess I should use EPSG:900913 as the base projection.


texts: How to filter by several parameters in GeoServer from the PostGIS table

texts:
On a layer from the GeoServer I linked a table to PostGIS, in that table I have a geometry column and several columns of related geometry information after which it should be filtered.
Now it should be when I request only certain parameters from the data bar, only that to bring back not the whole table.
I don't know how to pass those parameters to GeoServer and then filter them for me.
The code I've done so far:
var terrain = new ol.layer.Tile ({
          source: new ol.source.TileWMS ({
                    URL: "http: //127.0.0.1: 10000 / geoserver / ADMINISTRATION / WMS '
                    params: {'LAYERS': 'ADMINISTRATION: pdf', 'STYLES': 'LAND'},
                    format: new ol.format.GeoJSON (),
                    ratio: 1,
                    serverType 'geoserver'

          })
      });

In the table for example there is an identifier that is called city hall and has values from 1 to infinity now I would like it to show me in the field layer only the geometries that are in city hall 1.


texts: In QField is filtering one field based on values in another no longer working?

texts:
I am currently using QField (offline) 1.7.8, QGIS 3.16.1 and QField Sync 3.4.0 . QGIS is connected to a PostGIS database. I have created a form with a "Species" field that is filtered based on the values chosen in the "Genus" field. When I fill out the form in QGIS everything works, but when I try to use the same fields in QField in offline mode there is nothing, both fields are blank and give me no genera or species to choose from. In the past I have done this and QField worked flawlessly.
I am not sure if something has changed in QField and I need to do something different, or if this a bug. Has anyone else had this issue with recent versions of QField/Sync, or any suggestions for a solution?
**Some additional info I have one parent table and 2 child tables and the filtering issue is in a child table. Even when I open the "offline" project in QGIS that has been packaged for QField everything still works. The unfiltered lookup fields in the child table behave normally in QField. There must be something on QField's end causing the issue.


texts: Creating polygon from points with other attributes in pgRouting

texts:
I have table, which looks like this:

This is the result of pgr_drivingdistance.
I would like to create a polygons for each mz.

I tried to use group by.

drop table if exists service_area_polygon;
create table service_area_polygon as (
SELECT ST_SetSRID(pgr_pointsAsPolygon('SELECT id_node::integer as id, mz
, ST_X(the_geom)::float AS x, ST_Y(the_geom)::float AS y
FROM service_area_points group by mz, id_node, the_geom'),2180));

But this query creates one polygon with a extent of all points.
Something like this:

I also tried to use partition, but I'm new to these issues.

The new table with polygons should have the mz attribute.

I hope, that my problem is clear.


texts: FME workspace runner mangles WFS parameters

texts:
When running a workspace using a WorkspaceRunner transformer in FME some non-standard query string parameters in the WFS url are discarded and/or seemingly randomly replaced by equals signs.

Setup is like this:

Workspace "control":
  creator -> workspace runner (this runs the workspace "wfsread", parameter WFS URL reads https://[serverpath]?SERVICE=WFS&REQUEST=GetCapabilities&NodeId=287

Workspace "wfsread":
  wfs reader -> shapefile writer
  reader URL is "https://[serverpath]?SERVICE=WFS&REQUEST=GetCapabilities&NodeId=287"

Now, when running the control workspace, the log of workspace "wfsread" starts with

  2019-01-15 10:12:42|   0.0|  0.0|INFORM|FME Configuration: Command line arguments are 'C:\FME\fme.exe' 'D:/wfsread.fmw' '--SourceDataset_WFS' 'https://[serverpath]?SERVICE=WFS=GetCapabilities=287' '--FEATURE_TYPES' 'L525 L50 L48 L383 L382 L381 L377 L34' '--DestDataset_ESRISHAPE' 'D:\'

...so FME is somehow throwing away the characters "NodeId" and replacing ampersands with equal signs in the call to the wfsread workspace.

Is this a known issue, and/or are there workarounds/settings that can mitigate this?


texts: Alternate of Near Distance Tool of ArcGIS 10.1

texts: I have basic license level of ArcGIS 10.1, and near distance tool is not available with it. Can anyone please suggest of another similar purpose tool, under basic license, to get distance of one polygon from another polygon or any workaround?


texts: Feature service is not added to map using ArcGIS API for JS?

texts:
I'm trying to search features using the "esri/search" module. I'm passing the parameters to callback function (Search) to add a feature layer and search features, but i don't see any features added on map in first place. Please check below the function call:

var mysources = mySearch.get("mysources");
mySearch.on("load", function(){
  mySearch.Set("mysources",mysources);
  mysources.push({
    myfeatureLayer: new FeatureLayer("http://arcgis-serveruf-1981283418.us-east-1.elb.amazonaws.com/arcgis/rest/services/MyFiles/LandParcels/FeatureServer/0"),
      mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
      searchFields: ["ADDRESS"],
      suggestionTemplate: "Address : ${ADDRESS}",
      exactMatch: false,
      outFields: ["*"],
      name: "land parcels",
      labelSymbol: textSymbol,
      placeholder: "address of parcels",
      maxResults: 6,
      maxSuggestions: 6,
      enableSuggestions: true,
      minCharacters: 1
  });
  map.addLayer(myfeatureLayer);
});


texts: Why is Oracle Spatial Query Layer Not Displaying in ArcMap?

texts:
I am having an issue with an Oracle query layer using SDO_Geometry not displaying.

I have created a polygon layer/table within oracle, registered the table with the SDO_GEOM_METADATA and created a spatial index...twice.

Unfortunately, I am still unable to see the polygons display in an ArcMap query layer. I am able to bring the data into ArcMap, open the table, select and identify it and actually see the shape flash, so I know it is there.

Also, I can export a shapefile and the shapes render perfectly fine, but for some reason they will not display from the query layer directly.

I have checked over the index parameters in the registry and nothing seems off, see:

INSERT INTO all_sdo_geom_metadata
            (owner,
             table_name,
             column_name,
             diminfo,
             srid)
     VALUES ('DB',
             'TABLE',
             'GEOM',
             MDSYS.sdo_dim_array (MDSYS.sdo_dim_element ('X',
            -180,180,5.0e-10), MDSYS.sdo_dim_element ('Y',-90,90,5.0e-10)),4326);
index creation:
CREATE  INDEX wm_dim_table_si ON TABLE (geom) INDEXTYPE IS MDSYS.spatial_index
PARAMETERS ('tablespace=DB_INDEX03  layer_gtype=MULTIPOLYGON');

I've done everything, others have prescribed but I still can't view the data directly...
Any hints as to what else I should try or should be looking at? Any help is greatly appreciated.

My nexts thoughts are that I should check and validate the geometry on the polygons.
Also, there are about 30,000 polygons in this table, is there a threshold in ArcMap for query layers? It wouldn't seem so, however I am at a loss right now...


texts: Replacing non-English characters in attribute tables using ArcPy and Python?

texts:
I have a few shapefiles where some of the attributes contain the non-English characters ÅÄÖ. Since some queries doesn't work with these characters (specifically ChangeDetector), I tried to change them in advance with a simple script and add the new strings to another field.

However, change in characters works fine but not updating the field with arcpy.UpdateCursor.

What is an appropriate way of solving this?

I have also tried to do this via the Field Calculator while posting "code" to the codeblock, with the same error.

Error message:
Runtime error
Traceback (most recent call last):
  File "", line 1, in
  File "c:/gis/python/teststring.py", line 28, in
    val = code(str(prow.Typkod))
UnicodeEncodeError: 'ascii' codec can't encode character u'\xc4' in position 3: ordinal not in range(128)

Code:

# -*- coding: cp1252 -*-
def code(infield):
    data = ''
    for i in infield:
##        print i
        if i == 'Ä':
            data = data + 'AE'
        elif i == 'ä':
            data = data + 'ae'
        elif i == 'Å':
            data = data + 'AA'
        elif i == 'å':
            data = data + 'aa'
        elif i == 'Ö':
            data = data + 'OE'
        elif i == 'ö':
            data = data + 'oe'
        else:
            data = data + i
    return data

shp = r'O:\XXX\250000\DB\ArcView\shape.shp'

prows = arcpy.UpdateCursor(shp)

for prow in prows:
    val = code(unicode(str(prow.Typkod), "utf-8"))
    prow.Typkod_U = val
    print val
    prows.updateRow(prow)

The values of Typkod are of the type:
[D, D, S, DDRÄ, TRÄ] etc.

I use ArcMap Basic (10.1) on Windows 7.

New Error message:
Runtime error
Traceback (most recent call last):
  File "", line 1, in
  File "c:/gis/python/teststring.py", line 29, in
    val = code(unicode(str(row.Typkod), "utf-8"))
UnicodeEncodeError: 'ascii' codec can't encode character u'\xc4' in position 3: ordinal not in range(128)

>>> val
'DDRÄ'
>>> type(val)
type 'str'

It appears as if the output from the function is wrong somehow. When there's ÅÄÖ involved it returns data = u'DDR\xc4' and not (as was my intention) data = 'DDRAE'. Any suggestions on what might cause this?


texts: I have a set of lat and long coordinates whihc i have uploaded to be a table in mapinfo (the latest one) as you can tell i'm not much of a GIS wizz. I can't seemingly plot these coordinates in London they always plot in the sea somewhere below the ivory coast - this means my longitude is correct but my latitude isn't. However, coordinates exactly the same had been plotted by my colleague before,on the same layer successfully in london. Is there something i'm missing? Here is an example coordinate: 51.475566    -0.201027


texts:
I have a DEM and I want to fill every sink with a Value of 2 meter. If a sink is less than 2 meter, it is filled full and if it is bigger, it is filled until it reach 2 meter.

The problem is that this does not work with the Spatial Analyst Tool Fill. If the sinks are deeper than the Z-Limit, they will not be filled.

After the fill, the original DEM is subtracted from the filled DEM.
With this, I want to get a sink raster, which will show the depth of the water if the sinks are filled up to 2 meters.

I also tried the answer of the question asked here (Partially fill a sink - ArcGIS), but it does not really work.

This image shows a sink of my dem. The sink is roughly 10 m deep.

Here you can see the area after the method used in (Partially fill a sink - ArcGIS). The dark blue colour has the value 2 and the white 0. That cannot be, because the outer to inner part of the sink has a difference of 10 meters.

This image shows how the raster roughly should look like (I am sorry, I´ve got only Paint at my disposal). You can see clearly how deep the water is at each spot and is filled up to 2 meters.


texts: Clipping Imagery Services Layers from ArcGIS Online

texts: The new Landscape layers such as the various Soil type imagery services look like a great resource. Does anyone know if it possible to clip out small areas of the Landscape imagery service since I only need the soil data on a small area.  Clipping it will look much cleaner if I can clip to the feature class of my area of interest. When trying to run the Clip tool in ArcMap 10.2, I can't get the imagery service to show as an option to clip.


texts: SQL 'insert into table(values)' using QGIS-modeller

texts:
How can I run a precalculated SQL query using QGIS-modeller to insert values in an existing PostgreSQL table?

I know I can use import layer in DB-manager to load features in an existing table, but I'm working on a more complex model to import features in different related tables.

For testing purposes I created a simple table called test with 2 fields (name and type) and an autogenerated id-field:

I have a .dbf with the same attributes that I want insert to the DB-table:

I made a model that:

Refactored the attributes into one field with following expression:

concat(' (','\'', "name", '\',',"type", ')')

Which gives me this result:

Execute SQL to combine all the features into one:

select 'insert into test (name, type) values ' || group_concat(sql) || ';' as sql from input1

with result:

PostgreSQL execute SQL to run the calculated SQL query, using a pre-calculated value attribute(@Execute_SQL_OUTPUT, 'sql' ) (because algorithm output does not work)

This gives me an error saying:

  can't execute an empty query

Running the query 'by hand' within DB-manager works fine.  So the query is not the problem.

Adding the query into the 'PostgreSQL execute SQL-algorithm' as a value does the job.

So, how can I get the result from the previous query into 'PostgreSQL execute SQL-algorithm'? Or is it a bug?

more details:

QGIS version: 3.10.2-A Coruña
QGIS code revision: d4cd3cfe5a
Qt version: 5.11.2
GDAL version: 3.0.3
GEOS version: 3.8.0-CAPI-1.13.1
PROJ version: Rel. 6.3.0, January 1st, 2020
Processing algorithm…
Algorithm 'insert_into_db' starting…
Input parameters:
{ 'input' : 'C:/Users/Pieter/Desktop/GIS/insert_into_db/input.dbf|layername=input', 'qgis:executesql_1:sql_output_concat' : 'TEMPORARY_OUTPUT', 'qgis:refactorfields_1:output1' : 'TEMPORARY_OUTPUT' }

Prepare algorithm: qgis:refactorfields_1
Running Refactor fields [1/3]
Input Parameters:
{ FIELDS_MAPPING: [{'expression': 'concat(\r\' (\',\r\n\'\\\'\', "name", \'\\\',\',\r\n"type", \r\n\')\'\r\n)', 'length': 250, 'name': 'sql', 'precision': 0, 'type': 10}], INPUT: 'C:/Users/Pieter/Desktop/GIS/insert_into_db/input.dbf|layername=input', OUTPUT: 'TEMPORARY_OUTPUT' }
OK. Execution took 0.012 s (1 outputs).
Prepare algorithm: qgis:executesql_1
Running Execute SQL [2/3]
Input Parameters:
{ INPUT_DATASOURCES: ['NoGeometry?field=sql:string(250,0)&uid={d628877a-93e0-4c53-80d8-39fc2d902fd1}'], INPUT_GEOMETRY_CRS: None, INPUT_GEOMETRY_FIELD: '', INPUT_GEOMETRY_TYPE: 1, INPUT_QUERY: 'select \'insert into test (name, type) values \' || group_concat(sql) || \';\' as sql from input1', INPUT_UID_FIELD: '', OUTPUT: 'TEMPORARY_OUTPUT' }
OK. Execution took 0.234 s (1 outputs).
Prepare algorithm: qgis:postgisexecutesql_1
Running PostgreSQL execute SQL [3/3]
Input Parameters:
{ DATABASE: 'postgis_25_sample', SQL: None }
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\tools\postgis.py", line 840, in _exec_sql
cursor.execute(sql)
psycopg2.ProgrammingError: can't execute an empty query

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\algs\qgis\PostGISExecuteSQL.py", line 71, in processAlgorithm
db._exec_sql_and_commit(str(sql))
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\tools\postgis.py", line 852, in _exec_sql_and_commit
self._exec_sql(c, sql)
File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\tools\postgis.py", line 843, in _exec_sql
+ e.cursor.query.decode(e.cursor.connection.encoding))
AttributeError: 'NoneType' object has no attribute 'decode'

Error encountered while running PostgreSQL execute SQL
Error encountered while running PostgreSQL execute SQL
Execution failed after 0.47 seconds

Loading resulting layers
Algorithm 'insert_into_db' finished


texts:
I'm trying to find records from two tables with matching point geometry. The tables have different EPSG types. I do this:

SELECT count(*) FROM t1
WHERE st_intersects(t1.geom, (SELECT st_transform(t2.geom, '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', 3008) FROM t2));

but get

  ERROR:  more than one row returned by a subquery used as an expression

I've also tried:

SELECT count(*) FROM t1
WHERE st_intersects(t1.geom, (SELECT st_transform(t2.geom, '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', 3008) FROM t2 LIMIT 1));

and get 0 rows.

I'm unsure now if there are no matches (it is a possibility here) or if my SQL requries improvement?


texts: Where to see smallest size of segmented object in eCognition?

texts: I used ecognition to conduct segmentation, and I'm wondering where I can see the segmentation result, like the smallest size of object segmented.


texts: Merge an entire Geodatabase in ESRI ModelBuilder

texts: How can I merge an entire Geodatabase in the ESRI (ArcView) ModelBuilder? I have several FeatureClasses created into the GDB and in the model I plan to merge them all. Right now I have imported all the FCs by hand into the Model. This is not the solution I am aiming for because the names and the amount of FCs could vary.


texts:
I am trying to detect cloud shadows in Landsat images.
I came up with producing good shadows and cloud masks. Now I'd like to implement my cloud mask adding cloud shadows.

Do You have a good idea to do that?

OPTION 1: I thought combining the two masks I have to be a good idea, but I would need an oriented buffer to do that, given that I know the solar azimuth. I found a How to create an oriented buffer using arcpy? with a super good explanation, but I have no idea how to translate the R prototype code provided there into Python (which I use), as I have too few programming skills.

Does anybody know how to deal with this (maybe @whuber himself who provided the code)?

OPTION 2: As I know sun azimuth and angle (from .MTL file of Landsat), I thought to be a good idea to create a fake DEM with the cloud mask inside. Clouds would have their elevations as value in the DEM. Then I would cast the shadows against this "fake DEM" to obtain my cloud shadow mask (or to "play with" it and my shadow mask). Anyway, I don't know a good way to estimate cloud elevations, and, having tried with 5000, 8000, and 10000 meters values with very poor results, I am not totally convinced this is the best method.


texts: Extracting vertices in line segments using ArcGIS Desktop?

texts:
I have ArcGIS Desktop 10.6 and have many polyline features and want to be able to extract the vertices of each line segment that looks like this:

I want to keep the line ID associated and be able to link the points back to the polyline and also count and order the number of vertices along the polyline...

Ultimately I want to convert those points to polygons somehow select the inside polygons and delete them to create something like this:

I see other python scripts here but nothing that does what I need it to do. I've also tried feature to polygon but these FC's are exported out of a third party app and for some reason that tool doesn't work.

So what I have a polyline feature class that looks like a polygon but not closed. I want to extract the vertices along with assocaited line ID and number the vertices per line in order. Take the points shp and create a polygon.

Okay I created a model that works for what I needed...


texts: Rasterise shapefile as 1 band (grayscale)

texts:
My code converts shapefiles to rasters.
I had first done this in colour, but now I need greyscale (1 band).
The code seems to create this for the first raster file, but the other raster files seem to have no bands (previously they all had 3).

shapefiles = [folder + '/' + file for file in os.listdir(folder) if 'shp' in file]
output = 'output'
if not os.path.exists(output):
    os.mkdir(output)
RefImage = ref_image
gdalformat = 'GTiff'
datatype = gdal.GDT_Byte
burnVal = 1 #value for the output image pixels

# Get projection info from reference image
RasImage = gdal.Open(RefImage, gdal.GA_ReadOnly)

for Shapefile in shapefiles:
    ds = ogr.Open(Shapefile)
    ds_layer = ds.GetLayer()
# convert a shapefile to raster
def convertShpToRaster(inputfile, band, size):
    shp = ogr.Open(inputfile)
    lyr = shp.GetLayer(0)
    lyrname = lyr.GetName()
    #color = plt.gray()

    subprocess.call('gdal_rasterize -ts {0} {1} -l {2} -burn {3} ./{4} ./{5}.tif' \
                .format(size[0], size[1], lyrname, band, inputfile, output + '/' + lyrname))

def inputToShapefile():
    # input color, resolution and shapefile location
    #color = plt.gray()
    band = RasImage.GetRasterBand(1)
    size = RasImage.RasterXSize, RasImage.RasterYSize

    # call the function
    for inputfile in shapefiles:
        convertShpToRaster(inputfile, band, size)

inputToShapefile()

Can somebody advise? Should I rasterise as before and then drop bands instead?

I am using Python 3.7.


texts:
I got two shapefile of the same coordinate reference system (CGRS93_LTM), one in QGIS and the other in ArcMap. However they got different prj. files
QGIS file:

PROJCS["Transverse_Mercator",
 GEOGCS["GCS_WGS_1984",
  DATUM["D_unknown",
  SPHEROID["WGS84",6378137,298.257223563]],
 PRIMEM["Greenwich",0],
 UNIT["Degree",0.017453292519943295]],
PROJECTION["Transverse_Mercator"],
 PARAMETER["latitude_of_origin",0],
 PARAMETER["central_meridian",33],
 PARAMETER["scale_factor",0.99995],
 PARAMETER["false_easting",200000],
 PARAMETER["false_northing",-3500000],
 UNIT["Meter",1]]

ArcMap file:

PROJCS["CGRS_1993_LTM",
 GEOGCS["GCS_CGRS_1993",
  DATUM["D_Cyprus_Geodetic_Reference_System_1993",
  SPHEROID["WGS_1984",6378137.0,298.257223563]],
 PRIMEM["Greenwich",0.0],
 UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
 PARAMETER["False_Easting",200000.0],
 PARAMETER["False_Northing",-3500000.0],
 PARAMETER["Central_Meridian",33.0],
 PARAMETER["Scale_Factor",0.99995],
 PARAMETER["Latitude_Of_Origin",0.0],
 UNIT["Meter",1.0]]

How can i define Datum, ProjCS and GeogCS in QGIS in order to be the same?


texts: Create elevation-service from geotiff using python and gdal

texts:
I'd like to create a little server for us, that tells us the elevation of a point in a geotiff file. I already got it working basically using the following:

dataset = gdal.Open('test_file.tif', gdal.GA_ReadOnly)
geotransform = dataset.GetGeoTransform()
px, py = convertLatLngToPixels(geotransform, lat, lng)
val = self.dataset.GetRasterBand(1).ReadAsArray(px, py, 1, 1)

Now the issues is, that our real raster file, is a vrt that combines several geotifs. Using my code now on that vrt does not work any more...

Does anybody know any solution to using a vrt file?


texts: How to mix coordinate systems in GeoDjango?

texts:
This is a beginner question about working with different coordinate systems in GeoDjango.

I've got a GeoDjango object with a polygon field, stored as OSGB coordinates:

class Region(models.Model):
    name = models.CharField(max_length=30)
    geom = models.PolygonField(srid=27700)
    objects = models.GeoManager()

I'd like to be able to query it using standard WGS84 coordinates:

def api(request):
    lat = float(request.GET.get('lat'))
    lng = float(request.GET.get('lng'))
    pnt = Point(lat, lng, srid=4326)
    a = SettlementLocalRegions.objects.get(geom__contains=pnt)

But this gives me no results, even though I know at least one object should contain the point.

I figure I'm doing something wrong with my coordinate systems. Do I need to convert pnt to OSGB somehow? I thought the point of GeoDjango was that you could mix and match coordinate systems and it coped.


texts: Panning to specific feature using ArcPy?

texts:
I am attempting to use python and arcpy to automate the creation of some standardized maps I need to export for feature classes from various projects.  I have a set of points and I need to export a single map for each point in the feature class, showing only one point at a time.

I can iterate through the features and export the maps just fine.  However, the maps that come out are not panned to the right area, most of the time they are totally off.

I am using lyr.definitionQuery to hide all the features but one, then I need to pan my data frame so that it is centered on that feature, then export the map, and repeat for all the features in the feature class.  The best I could find to recenter the data frame was the df.panToExtent function.  Which works fine when I run it from the python window in ArcMap but doesn't seem to keep up when it runs inside the loop.  Any help is MUCH appreciated!

Here is my loop:

while isoPosition < isoCount: #iterate through the feature class
    query = "[Type] = 'Isolate' AND " + "[FieldID] = " + "'" + FieldIDs[isoPosition] + "'"

    #center on feature
    df.panToExtent(lyr.getSelectedExtent())

    #Refresh the view to load new image tiles
    arcpy.RefreshActiveView()

    #Set the layers definition query equal the query text
    lyr.definitionQuery = query

    #String for a unique file name for each map, corresponding to the current Field ID in FieldIDs
    filename = FieldIDs[isoPosition]

    #Export the map document:
    arcpy.mapping.ExportToPNG(mxd,filename)

    #Move forward in the FieldIDs array, so that next time through we export the next map
    isoPosition += 1


texts: Linear referencing to calculate maximum distance between points by group along line (river system) using ArcGIS Pro

texts:
I am using ArcGIS Pro.

I have a series of points representing animal locations that were tracked in a small stream. I would like to get the distance between the most upstream and downstream location for each individual (also known as the Linear Home Range), following the river path. Essentially the largest distance along this line between points, grouped by individual. I have downloaded a polylineZM file for the stream - and I understand that I can use Linear Referencing to get the distances between these points, and a measure of the point location in reference to the beginning of the line. But, I am looking for a little more than this, as this result gives me an output that I would need to manually manipulate to get my answer. If there is a code or workflow that would help me to easily calculate this without having to split the data set into a file for each individual, manually rerun it for each individual, sort output by the distance along the line and then manually subtract the largest and smallest points, that would be ideal.

The most recent point may not necessarily be the most upstream point, as the animals are not only traveling in one direction. Below is a screen grab of my points -

If there is a need to write an ArcPy script or a suggestion for manipulating the data after the Linear Referencing output file after it is created - I am all ears.


texts: Seeking road network of Sweden with speed limits

texts:
I need road network data of Sweden with the speed limit of each segment of the road. I downloaded it from https://download.geofabrik.de/ but unfortunately, there are a lot of missing values there.
Is there any other website to get this data from free?


texts: Transparent png tiles showing up weirdly in QGIS?

texts:
I have Ubuntu 14.04 running QGIS 2.6. I'm loading a raster layer, but the transparent PNG shows up poorly. The black lines along roads aren't showing up and the anti-aliasing looks off. I'm thinking that QGIS might be clipping the black values. I've fiddled with the raster settings but I can't get it to work. Here's the original file.

With resampling changed to Bilinear


texts: Installing GDAL with MrSID fails on Ubuntu

texts:
I'm trying to install GDAL from source on Ubuntu or CentoOS both, but it fails both times. I also tried GDAL 3.0 which did the same thing.

Is there any way easily get GDAL working with MrSID? I can install GDAL with apt install gdal-bin, but there is no MrSID support.

I'm using

./configure --with-mrsid=/home/centos/gdal-2.4.3/mrsid/Raster_DSDK
make

The output of make on Ubuntu is:

mrsidstream.h:64:17: warning: âvirtual const char* LTIVSIStream::getID() constâ can be marked override [-Wsuggest-override]
     const char* getID() const;
                 ^~~~~
mrsiddataset.cpp:214:23: warning: âvirtual LT_STATUS MrSIDProgress::setProgressStatus(float)â can be marked override [-Wsuggest-override]
     virtual LT_STATUS setProgressStatus(float fraction)
                       ^~~~~~~~~~~~~~~~~
../../GDALmake.opt:670: recipe for target '../o/mrsiddataset.lo' failed
make[2]: *** [../o/mrsiddataset.lo] Error 1
make[2]: Leaving directory '/home/ubuntu/gdal/gdal-2.4.3/frmts/mrsid'
GNUmakefile:15: recipe for target 'mrsid-install-obj' failed
make[1]: *** [mrsid-install-obj] Error 2
make[1]: Leaving directory '/home/ubuntu/gdal/gdal-2.4.3/frmts'
GNUmakefile:103: recipe for target 'frmts-target' failed
make: *** [frmts-target] Error 2

The output of make on CentoOS 7 is:

/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTIOFileStream::~LTIOFileStream()'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTIOFileStream::initialize(char const*, char const*)'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTIMetadataDatabase::~LTIMetadataDatabase()'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTIImageStage::read(LizardTech::LTIScene const&, LizardTech::LTISceneBuffer&)'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::MrSIDImageReader::create()'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTINavigator::LTINavigator(LizardTech::LTIImage const&)'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTINavigator::setSceneAsULWH(double, double, double, double, double)'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTIGeoCoord::get(double&, double&, double&, double&, double&, double&) const'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::LTISceneBuffer::LTISceneBuffer(LizardTech::LTIPixel const&, unsigned int, unsigned int, void**, unsigned int)'
/home/centos/gdal-2.4.3/.libs/libgdal.so: undefined reference to `LizardTech::J2KImageReader::initialize(LizardTech::LTIOStreamInf&, bool, float, bool)'
collect2: error: ld returned 1 exit status
make[1]: *** [gdalserver] Error 1
make[1]: Leaving directory `/home/centos/gdal-2.4.3/apps'
make: *** [apps-target] Error 2

I just downloaded the SDK from https://www.extensis.com/support/developers
I tried 5.3.1 and 4.8.2.


texts:
How do you get the cell size of a GRASS raster in QGIS?
QGIS has properties, but no raster cell size.


texts: Create GPS (gpx) routes from GeoPandas geodataframe

texts:
I have a GeoPandas geodataframe with LineString geometries and I'd like to export them into .gpx files as GPS routes.

I've searched for hints in the web but found anything useful.

Are there GeoPandas specific way to do this or from other packages like Shapely?


texts: Filter date from featurecollection based on geometry points on Google Earth Engine

texts:
I am creating a code to automate the process of downloading point-based reflectance data tables that comes from a JSON file. The code is almost ready and I can reduce and export a table containing the reflectance data per band at each point. The problem is that I'm generating a table with all the reflectance values for each point I've created for all available dates in the time interval between points. I need to generate a table with only data referring to the date the point was created.

Any idea how to do this?

This is the code working so far: Code here

var ict = ee.Dictionary({RAWJSONDATA})
var feats = ee.FeatureCollection(ee.List(dict.get('pings')).map(function(feat){
  feat = ee.Dictionary(feat);
  var geom = ee.Geometry.Point([ee.Number.parse(feat.get('longitude')), ee.Number.parse(feat.get('latitude'))]);
  var getDate = ee.String(feat.get('datetime'));
  var systemDate = ee.Date.parse("dd MMM YYYY hh:mm:ss aa", getDate).millis();
  return ee.Feature(geom, ee.Dictionary(feat)).set('system:time_start', systemDate);
}));
feats = ee.FeatureCollection(feats.setMulti(dict));

print(feats);

var rrs = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
    .filterDate('2018-2-22','2018-2-24') //HOW DO I FILTER WITH THE "systemDate" ??

    .select('Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555');

var reduceR = function(image) {
  var reducedVal = image.reduceRegions({
    collection: feats,
    reducer: ee.Reducer.first(),
  });
  return reducedVal;
};

var mRrs = rrs.map(reduceR);
mRrs = ee.FeatureCollection(mRrs).flatten();

Export.table.toDrive({
   collection: mRrs,
   description: 'aaaagoraaa',
   fileNamePrefix: 'vaiii',
   fileFormat: 'CSV',
   selectors: ('Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555'),
 });


texts: How to set a vector point on an exact lat/long coordinate

texts: I have a vector polygon I'm trying to edit, I'd like to use exact lat/long coordinates for one of the vector points. How can I do this?


texts:
I have 423 spot.vgt ndvi from first April 1998 till 31 Dec 2009.
I want to know how to stack these NDVI to have one image for time series analysis?
Or better to ask how is it possible to analyze a pixle change during these period?


texts: ArcGIS API for JavaScript won't add graphics to map

texts:
var point = new Point(-49.9, 39.3)
var symbol = new PictureMarkerSymbol('http://www.prww-weather.com/hurricane/icons/post-tropical-cyclone.png', 50, 50)
var graphic = new Graphic(point, symbol);
map.on("load", map.graphics.add(graphic));

These four lines of code use the ArcGIS API for JavaScript to add a point to the map. However, for some reason, the code doesn't work and I get this error in the console:

  Uncaught TypeError: Cannot read property 'add' of null

Can anyone tell me why this is occurring and how to fix it?


texts:
I am currently trying to draw a point with the ArcGIS Javascript API.

Pretty much like this example, http://help.arcgis.com/en/webapi/javascript/arcgis/demos/toolbar/toolbar_draw.html, except the point will not be drawn with the mouse, but from an javascript object.

If we look at the code, we see that a geometry object/array is passed in the function.

function addToMap(geometry) {
    toolbar.deactivate();
    map.showZoomSlider();

    var symbol = new esri.symbol.SimpleMarkerSymbol(esri.symbol.SimpleMarkerSymbol.STYLE_SQUARE, 10, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([255,0,0]), 1), new dojo.Color([0,255,0,0.25]));

    var graphic = new esri.Graphic(geometry, symbol);
    map.graphics.add(graphic);
  }

If someone could write the structure for that object, and a method to draw the point (such as a button, with an onclick event, containing fixed xy coords)

Dumping geometry looks like this http://pastebin.com/P9wfL36b (ouch!)

This code will be used in a GPS tracking experiement, where I'll pull the current position from a database every 20 seconds or so. (details will not be covered in this post)

Bonus question: How do I convert long/lat to xy coords from the sample?

I use UTM zone 29N, wkid: 32629


texts: Inserting/Update a polygon into a field with null geometry in Oracle

texts:
How can I copy a previously created polygon that is in a temporary layer, to a field with null geometry in Oracle?
I can view the field in the table, use the "Add Part" tool to create a new polygon, but I can't copy a previously created polygon. I use QGIS 3.14.1-Pi


texts: Seeking tool to use in ArcPy that does opposite of Near?

texts:
I have a shapefile with many (1000s) points grouped by attributes. I need to determine the distance between each point and the point furthest from it that shares the same group attribute.

Is there a tool that does the opposite of Near (proximity toolset) i.e. that would ideally operate like this Near By Group tool?


texts: Story Map SpyGlass without ArcGIS Online?

texts:
I am looking to see if there is an implementation example of the Spy Glass widget without ArcGIS Online or Story Maps.  I have found the LayerSwipe documentation and can implement that fine, but I have not found a way to do the Spy Glass version in my own implementation.

I am using it in a custom web app that will be hosted on our server and want to avoid using AGOL/ESRI styled items.  However, the Spy Glass caught someone's eye so I have to figure out how to implement it in our setup.

I am assuming it has to deal with using the LayerSwipe widget and different CSS styling/event handling (e.g. mouse move), but I am hoping to see an implementation so I don't have to do it all from scratch.

I have downloaded the source code from Github for the StoryMap template, but as I go through it, it seems like there is a lot more there than what I want.  I am searching through trying to pick and choose what code/CSS I want, but I am hoping someone else may have done something similar to save me some time.  My end result is going to be wrapping the widget in an AngularJS directive, so the simpler it is, the better.


texts: WFS GetCapabilities in human readable format?

texts:
I've got that impressions that GetCapabilites is essential part of machine to machine data interchange.

However, I also feel that the returned XML includes information that should be presented in human redable format or what?
For example list of available return formats which is quite crucial to developer to know what formats can requested?


texts: Cutting Raster using two Polygons as Masks

texts:
I am trying to cut this raster in order to drop the outside of the walls white area and the area with the dots.
I made one step using as mask the  green polygon but i couldn't remove the inside of the walls area with the dots.
Any idea?


texts:
Us lucky fellows in the UK have been able to access public geodata from Ordance Survey under an open licence since April this year.  Some of the datasets, including VectorMap, are broken into 100 files per grid square, each containing multiple layers, in this case in shape format.  To cover the wider London area with VM you need to load 400 x 12 or so layers, with each layer identically named.  Any tips on the best ways to do this, preferably using free or open source tools?

http://www.ordnancesurvey.co.uk/oswebsite/opendata/