texts: How to calculate the length of the landcover types a path crosses in QGIS

texts:
I have a path which represents a proposed pipeline route.  I need to find the length of the various landcovers the route crosses.  Not area within a buffer.  I found a ArcGIS help topic that kinda shows what I want to do: http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Overlay_analysis

Is this the best way to do it? ie- convert raster to vector, overlay vector and line?

EDIT- the land cover file is a raster file.


texts: Converting Raster datasets into layers using ArcPy/ModelBuilder?

texts:
I am new to ArcGIS and ArcPy (just started using them on Wednesday).

I am using a model to create a Point Density plot from a Shape File, with the output being a Raster dataset.

I want to convert the Raster dataset into a new layer and modify the properties.

Is there a function that converts Raster datasets into layers, or is there another way to do it?

It is clear the the model is just a starting point, and that I have to add Python code for it to do what I want it to do.


texts:
I'm just starting with GIS, and I'm looking for a routing service that will allow me to avoid given regions on a map, or to build my own router that has this functionality. The site OpenRouteService.org has this capability, but the back-end of this website seems to be down for the last few days. So I'm wondering, what would be your advice for doing this? Either a service that will plot from point A to B given the constraint to avoid a region, or how to build my own version. Or a working clone of openrouteservice? What I've tried is downloading the github page of openrouteservice here:https://github.com/GIScience/openrouteservice. The provided .html gives a working front-end, but I have to say I'm a bit  confused how to start getting the back-end working.

I should mention I'm looking for something that will work in the U.S. northeast. Openroute seems to not have the U.S. in its coverage area currently.


texts:
I need to calculate possible point location based on its speed and some other factors. For example, there is a man in the middle of the forest and I need to calculate his possible location at some time periods. There are few facts about his movement:

His speed differs at certain periods of time - it can be raining, or it can be super hot - speed lowers / increases;
Some objects in space - swamps, rocks - can lower his speed as well;
At some time periods the man does not move at all.

So, here is the question - is there any algorithm that can calculate possible location of the man based on following parameters? The output should be a polygon for every given time_period.

{
  time_periods: [{
    from: 2012-12-12T20:00:00,
    to: 2012-12-12T22:00:00,
    speed: 5
  }, {}],
  objects: [{
    geodata: POLYGON(...),
    speed_modifier: 0.5
  }, {}]
}

Generated polygons will be visualized on the map afterwards.


texts: Create a sript for "point sampling tool" of QGIS

texts:
Is there any way to create a python script to the "point sampling tool" plugin from QGIS?

I have a very large number of Netcdf files which I want extrat the values from pixels to a irregular points shapefile.

Any idea?


texts: Constant future warnings with new pyproj

texts:
When I work with different CRS in GeoPandas, pyproj module always fires many FutureWarnings.

/home/culebron/.local/lib/python3.7/site-packages/pyproj/crs.py:77:
FutureWarning: '+init=<authority>:<code>' syntax is deprecated.
'<authority>:<code>' is the preferred initialization method.

I see them up to 5 in a row.

Here's a code that causes it:

In [1]: crs = 'epsg:3857'
In [2]: import geopandas as gpd
In [3]: df = gpd.read_file('rodovias.gpkg')

In [4]: df.crs
Out[4]: {'init': 'epsg:4326'}

In [5]: df.to_crs(crs)
/home/culebron/.local/lib/python3.7/site-packages/pyproj/crs.py:77: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method.
  return _prepare_from_string(" ".join(pjargs))
Out[5]:
                     name                                           geometry
0                Arco Sul  MULTILINESTRING ((-5487151.115 -2940693.219, -...
1                   Leste  MULTILINESTRING ((-5482946.592 -2929361.487, -...

As you can see, I have no syntax like '+init=...' in my code.

Is there a workaround to silence the warnings, or ideally to fix their cause?

(edit: and please, don't vote close. I did Google before asking this, and found absolutely nothing on topic. Bothering the devs at geopandas github was a bad idea too.)


texts: How do I install GRASS to run in QGIS 2.0.1?  GRASS is not listed in the Plugin Manager and attempting to access it via the Processing Toolbox gives me the following error:


texts: r.flow in QGIS 2.16.1 generates wrong output

texts:
I try to run the r.flow module to get the flow path length of my layer out.
The input Layer is a DEM with pixelsize 0.0083 and the selected CRS (EPSG:32646, WGS 84/ UTM zone 46).
First of all I generated a depressionless elevation with saga:fill sinks.
If I start the function, it runs and I get the three Outputlayer, but they are wrong. -->  pictures see below

My Data: Win 10, Qgis 2.16.1 with Grass 7 Plugin.
I have the same problem if run r.flow in Grass.

Which additional steps I have to do ?


texts: Is it better to calculate area of a global raster using long/lat or using an equal area projection (e.g. Galls-Peters)?

texts:
I am trying to calculate the area of values in a global spanning raster with values of NA, 1, 2 or 3. I would like to know the total area of each category, however, the map is currently in lat long (unprojected) coordinates.
Is it better to calculate the area using the unprojected coordinates (e.g. with the area() function in the raster package in r) or is it better to first project the raster to a projected coordinate system (e.g. Gall's Peters) and then calculate the area?


texts: ST_MakeLine equivalent for CIRCULARSTRING in PostGIS

texts:
I'd like to make a CIRCULARSTRING in PostGIS 2.0.0 out of points that are calculated as part of the query (rather than constants). For a LINESTRING I can do this using ST_MakeLine(). Is there an equivalent for CIRCULARSTRING or some other way to do this?

I have 3 points that are calculated in some way (in the example let's just use ST_GeomFromText):

SELECT ST_GeomFromText('POINT(10 10)'), ST_GeomFromText('POINT(15 10)'), ST_GeomFromText('POINT(10 20)')

What I want to get is CIRCULARSTRING(10 10, 15 10, 10 20).

If I run:

ST_LineToCurve(ST_MakeLine(ARRAY[ST_GeomFromText('POINT(10 10)'), ST_GeomFromText('POINT(15 10)'), ST_GeomFromText('POINT(10 20)')]))

That gives me LINESTRING(10 10,15 10,10 20), not a CIRCULARSTRING.

Same question for CURVEPOLYGON - what can I use instead of ST_MakePolygon (which returns "Shell is not a line" if I pass in a CIRCULARSTRING)?


texts: Maintain color range in GeoTIFF when cropping

texts:
I have a series of scanned maps stored as GeoTIFF raster images.  When I crop them, the color representation changes, and it appears that QGIS calculates a minimum value for each band and defaults to a Contrast Enhancement of "Stretch to MinMax", which darkens the image and makes it less readable.  I can correct the image by either changing the Min values in each band to 0, or by changing the Contrast Enhancement to "No enhancement," but these changes are stored in the project rather than in the raster files.
Is there a way to save the files so the color enhancement is stored with the raster file rather than with the project?
In the image below, the lower half shows a layer with the settings displaying similar to the original, the upper half shows the default with the cropped image.


texts: displaying scale of dyamic map service layer using ArcGIS JS API

texts: In ArcGIS JS API, I would like to know and display the scale of a dynamic service layer by using the onExtentChange event, but the LOD argument is only valid for tile map services. Is there an equivalent way for dynamic map services?


texts:
I am using QGIS 2.0. I have downloaded USGS 1:100,000 topographic maps which show adjacent areas. The topographic maps were digitized with the borders (I've seen them called "collars", I think) as part of the rasters.
When I load the vertically (north-south) adjacent rasters, the border areas of course overlap the map areas of each other.
I have tried to convert one raster to a vector so that I could use the vectorized map as a mask to crop the other raster, but it didn't work. (My computer just kept spinning -- sort of froze up.)
I have not found, here or elsewhere, discussions of using one raster to crop or clip another.

If I can use one raster map to crop another, adjacent-but-overlapping raster map, then how can I do it?


texts: QGIS failing to complete clip or intersect

texts: I have been trying to perform a intersect (then tried using clip) on a rather large dataset on QGIS. I have been able to complete similar geoprocessing using slightly smaller datasets and the layer. For this layer no errors are being reported but the clip/intersect function seems to be unable to move past 28% complete (tried multiple times). Does anyone have any advice or help they can offer beyond trying to maximize my computers available processing capability?


texts: Arcpy - AttributeError: PageLayoutObject: Error in executing ExportToPNG

texts:
I’m building a custom tool, where you input the city and state and the tool outputs an image file of a map with certain features buffered, labeled, etc.  I’m running into problems on the last step, exporting to PNG. Here’s the relevant code:

#Set workspace
import arcpy
import arcpy.mapping
import os
arcpy.env.overwriteOutput = True

yourstate = arcpy.GetParameterAsText(0)
city = arcpy.GetParameterAsText(1)
distances = arcpy.GetParameterAsText(2)
folderinput = arcpy.GetParameterAsText(3)
folderinput2 = '"'+ folderinput + '\"'
arcpy.env.workspace = folderinput2

(I deleted a bunch of code here…it’s for editing the mxd based on what state/city you select. FYI the GetParameterAsText(3) input is something like
C:\Users\first.last\Desktop\ArcGis\tool)

#exportasapng
mxd = arcpy.mapping.MapDocument("CURRENT")
fullname = '"' + city + yourstate + '.png' + '"'
arcpy.mapping.ExportToPNG(mxd, fullname, resolution=200)

Any suggestions? Here is the error when I try to run the tool:

Traceback (most recent call last):
  File "C:\Users\first.last\Desktop\ArcGis\tool\toolv6.py", line 79, in <module>
    arcpy.mapping.ExportToPNG(mxd, fullname, resolution=200)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1228, in ExportToPNG
    layout.exportToPNG(*args)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\mixins.py", line 586, in exportToPNG
    return self._arc_object.exportToPNG(*args)
AttributeError: PageLayoutObject: Error in executing ExportToPNG

If I replace
arcpy.mapping.ExportToPNG(mxd, fullname, resolution=200)

with something like the following the tool exports successfully:

arcpy.mapping.ExportToPNG(mxd, "BuffaloNY.png", resolution=200)

….which is frustrating since “BuffaloNY.png” is what is being passed through as the parameter (it’s the result for print = fullname).

I’ve tried lots of minor changes, making sure the arcpy.mapping.ExportToPNG parameter is being passed as a literal string, changing the direction of the slashes, putting in a mxd.save before exporting, inserting the full pathname instead of relying on setting the environment, etc. and nothing is working. From searching around online, it seems like this is a problem other people have run into, but the solution wasn’t readily apparent.


texts:
I need an advice for a service.

On our website we have a search for events. Users are especially interested in the location of event so they use that filter a lot but the problem is that in our database we have sometimes England sometimes UK so we cannot provide them the information they are looking for.

Do you have any experience with this? We store the location of events in plain text (e.g. "London, UK") and we would like to use something more sophisticated. So in the future user will be able to search for events in London and in 100km radius for example. (so we show them also the events in Oxford.)

I have no experience with GEO location so maybe it's something simple maybe not.


texts: How to delete empty fields from a Shapefile with QGIS or ogr2ogr?

texts: I have a Shapefile with a lot of empty columns but many of them are empty. Using QGIS or ogr2ogr, how do I remove all attribute columns where everything is NULL?


texts: Problems using ATCOR for Erdas Imagine '13

texts:
An output from ATCOR processing (besides an atmospheric correction and dehazing) is expected to be an image containing true reflectance values (instead of raw input DNs). As far as I know such data are commonly stored as a floating point values (and not 8 or 16 bit data), because each pixel represents a fraction (from 0 to 1).
   But each time I process an image (LS8, QB, etc) I get an output of the same type as an input image was (8b or 16bit data). Does anyone know why is that? (I used a proper calibration file each time, and also there is no output data type options dialog - it just prompts for a filename).


texts: Getting Id of selected feature in QcomboBox of PyQGIS?

texts:
I have a QcomboBox that shows attributes from a field, each attribute occurs only one time. I am trying to get the id of the selected attribute\feature so I can then use setSelectedFeatures(), and then zoomToSelected().
I though there are similar examples, I can’t figure out how to get only the id of the selected feature.

Does anybody have a suggestion?

here is my code:

def run(self):
    """Run method that performs all the real work"""
    layer = iface.addVectorLayer("C:\Users\carmel.han\Documents\QGIS\qgis_sample_data\shapefiles/regions.shp",
                                 "regions", "ogr")
    if not layer:
        print "Layer failed to load!"

    Feature_list = []

    self.dlg.comboBox.clear()
    for feature in layer.getFeatures():
        name =feature ["name_2"]
        Feature_list.append(name)
    self.dlg.comboBox.addItems(Feature_list)

    id=
    layer.setSelectedFeatures(id)

    iface.mapCanvas().zoomToSelected()


texts: Styling GeoServer pbf vector tiles in Leaflet

texts:
I'm serving up vector tiles in pbf format in GeoServer 2.11.2 (using the MVT extension) and consuming them in Leaflet 1.2.0 with Leaflet.VectorGrid 1.2.0 but am having trouble styling them.
The URL format I am using is:

http://xxxx/geoserver/gwc/service/tms/1.0.0/GIS:LSOA@EPSG%3A900913@pbf/{z}/{x}/{-y}.pbf

The vectors appear and render nicely but with a default blue line style. The docs state you should set the style for each layer within the vector tile set like so:

var vectorTileOptions = {
vectorTileLayerStyles: {
    // A plain set of L.Path options.
    water: {
        weight: 0,
        fillColor: '#9bc2c4',
        fillOpacity: 1,
        fill: true
    }
};

i.e. in this example 'water' is the layer name.

I am unsure what my layer name is because whatever I have tried, the style remains the default blue.

In GeoServer my layer title is "LSOA" and name is "GIS:LSOA" so I have tried the following without success:

LSOA: {
        weight: 0,
        fillColor: '#9bc2c4',
        fillOpacity: 1,
        fill: true
    }

JS Error with the colon

GIS:LSOA: {
        weight: 0,
        fillColor: '#9bc2c4',
        fillOpacity: 1,
        fill: true
    }

Nope:

pbfLayerOpts["pbfLayerStyles"]["GIS:LSOA"] = {
        weight: 0,
        fillColor: '#9bc2c4',
        fillOpacity: 1,
        fill: true
    }

All attempts result in the same default style.
I have also tried creating a layer group in GeoServer and adding one or multiple layers but no luck there. Also tried adding a layer group without a workspace (the GIS prefix) but in that case GeoServer returned an error to the pbf request.

If I use the OpenMapTiles service then everything works as expected so it seems the issue is with GeoServer/the extension and how it names layers/serves the data.

Has anyone got this set-up working before?


texts: How to calculate areas for data that covers the whole of Australia

texts:
I have a polygon shapefile showing mining areas across all of Australia. I want to calculate the area for each of the polygons however they have been provided in a Geographic Coordinate System (GDA_94).

Can anyone tell me what projected coordinate system I could use that would enable me to calculate the areas of each polygon?

Will I need to transform the data and how might this affect the calculations given the size of Australia (i.e. it covers multiple GDA or UTM zones)?


texts:
I just started working with the OpenGeo Suite.  I'd really like to use it for a County web map, but have a couple questions...

1) Is there a way to add a lat/lon or other X/Y location to the URL and have the map center and zoom to the point?  Something like "?lat=46.265&lon=-90.206&zoom=7" at the end of the URL.

2) Is the Query tool available in a published map?  It's not in the Choose Tools window when publishing.

Thanks.

=== UPDATE ===

I will drop part 2 and stick with part 1.

I guess I forgot to mention the crucial bit that I'm working in GeoExplorer.  I did have the tag, but don't fault anybody for missing it.

The bounding box isn't working in GeoExplorer.  I have the map published in an iframe.  Here's a screenshot of the test page without the bbox... i.imgur.com/Wc022Ee.jpg.

Here's a screenshot with the bbox... i.imgur.com/7MES4cn.jpg.  That's just a blank map that hides behind the real map and only shows up when the real map doesn't.

I do have a correct bbox... i.imgur.com/MbAOgVD.jpg in the code.

(I guess I need a higher rep before I can add more direct links, so you'll need to copy paste the image URLs if you want to see them.)

Here's the source behind the page...

Iron County Test Map

This is a test map, version T-001.  A disclaimer might go here.

I tried putting the "?bbox=" at the end of the "src=.../1".

Thanks for any help.  And thanks for the welcome.


texts: How to merge many individual KML layers into one?

texts: I'd like to open a USGS KML file of bedrock geology for an entire county, but there are over 2,000 individual polygon layers.  Does anyone know if there is a way to merge all of these layers together into one file before opening in QGIS?


texts:
I've created my own data file using gps information and loads of sensor data.
I'd like to create a gps file which will add this information to the geographic location. Unfortunately most common gps file types don't support additional data like the gpx file format for example.

Sensor data would be telemetric data like humidity, temperature, etc.

So my question is:
What is best practice in order to create a unified gps file with the possibility to add additional data.

I've created my own data file using gps information and loads of sensor data. I'd like to create a gps file which will add this information to the geographic location. Unfortunately most common gps file types don't support additional data like the gpx file format for example.

Sensor data would be telemetric data like humidity, temperature, etc.

So my question is: What is best practice in order to create a unified gps file with the possibility to add additional data.

Additional Info:
Video length: 6:18:05 hrs (378:04,500 Seconds)

Data record first timestamp: 10:03:25 o'clock
Data record last timestamp: 15:16:03 o'clock
==> Data record length: 5:12:38 hrs (312:31.000 Seconds)

Here are the pont in times for three different sync points.

     1. Take-Off   =>   2.Explosion (balloon) =>    3.Landing
video   23:53:352         184:19.886                  241:20:765
data    15:50.105         149:15.014                  196:36.998
offset  -08:03.183        -35:04.872                  -44:43:767

As you can see, the offset is getting bigger the longer the video/data runs.


texts:
I have ArcGIS Server Basic Enterprise Version 10.3.1 and a Enterprise Geodatabase in Oracle 12c.

I want to create a service which will serve data from a view as tiles (PNG, JPG, whatever).

What type of service should I create?

Edit for more information:

I need to display a road network on top of a baselayer (google maps, OpenStreetMaps, etc) in a Web Application (.NET Core) using the ESRI Leaflet library. My netowrk consists of about 20000 links and it is stored in an Oracle database table.

I have already made a Feature Access Service, but because of the magnitude of data, the resulting response is huge (over 20MB) and cannot be handled by the browser adequately.

Since I have accomplished a similar task in the past using the Open Source MapServer, I am thinking about the following approach with ArcGIS Server:

Instead of returning a huge geoJson response, I think it would be better to configure the service in order to return tiles (PNG/JPG/whatever) with the links, therefore mitigating the load.

While the road network will always be the same, the "coloring" of the links will vary according to user preferences. For example, one user will might want to "color" the links according to "Vehicles Per Hour" during Summer of 2014, and some other user will might want to "color" it according a projection for Winter of 2020 (tha calculations will be performed by another module and the results will be stored in Oracle. I have already created a database View for properly quering the database. On this view I have based the existing Feature Access Service).

So, to cut a long story short, I cannot create a cached map beforehand; I'd rather like to create some sort of "dynamic" tiles, upon each request. This is something that I know (by first-hand experience) MapServer can do easily. Can ArcGIS Server Basic Enterprise behave the same way too?


texts: How to Set Page Numbers with Strip Map Index Features Tool

texts:
I have used the Strip Map Index Features tool to create an Index Layer for Data Driven Pages.  The problem is with the page numbering that the tool creates.  Some of the page numbers are in order, but others seem random.  Here is the tool as I ran it:

Here is the result I got:

How can I modify the tool so that the page numbers are in some order?  My preferred order would be along individual polylines from West to East and North to South.


texts: How to keep 0 decimals using QGIS Field Calculator's "to_string"

texts:
Why does Field Calculator's to_string function rounds decimal numbers if they have zeros?
For example: to_string(70.0) gives 70 (textual). But I would need it to be 70.0.
And as I have to convert an entire column which might have other values such as 70.1, adding the .0 manually is not going to be a solution to my problem.
Any idea?


texts: Is there a way to extract z values from a DEM into a polyline and/or polygon featurclass in arcmap 10? I would like to have 3d polylines and polygons to make cross sections in a CAD, but currently I have no way to project to a surface.


texts: How to create feature layer in Google Maps?

texts:
In ArcGIS we have the ArcGIS Online portal kind of application where we can create a feature layer.

Is there any application where I can maintain the feature data in Google?

Or how we can do this in the Google Maps API?


texts:
I'm creating regression lines using matlab from points having UTM geographic coordinates.

I would like to import those regression lines into a map of the same area in ArcMAP.

My question is what kind of data do I need to extract from Matlab to be read and imported as polylines in ArcMAP?

Thank you for your help


texts:
In QGIS 3.6.1 Noosa I have a column in the attribute table (Install Date) formatted as DATE that I want to add 100 or 60 years (Replacement Date) to, also in DATE format, for use in Time Manager.

I have tried creating a new field in the Field Calculator but get NULL values out. Expression used was todate("Install_Da"+'100 Years").

I can get a result using Integer formats using left(Install_Da,4) then another column adding 100 years but then it doesn't run in Time Manager due to wrong format.

What is the proper way to add years to a date?


texts: why is my SDE feature class not showing up in arcmap

texts: Just wondering if anyone has had a similar problem. We have a feature class on an SDE with about 1500 polygons. The attribute table says everything is there, but when we try to view it in arcmap, nothing shows up. If you run a function, clipping for example, it doesn't recognize any geometries from the feature class. We've tried running repair geometry, but that didn't work. Is it something to do with the extent?


texts: Bulk conversion and import of OS MasterMap .gz files to .shp in ArcGIS Desktop

texts:
I have 91 .gz files procured from the Ordnance Survey and would like to import them all to ArcGIS as shapefiles. I know that the Quick Import tool in the Data Interoperabilty extension performs the process correctly for each file (using these steps), but it would be too tedious to do that for the 91 gz files into 91 file geodatabases.

How can I automate this process so that the files are converted and each subsequent feature class is merged into the resultant layers in one file geodatabase?

I'm in my early days of python scripting.

Inside the .gz files, there are GML files. So the question is how to batch convert these GML files to a single feature class (e.g. inside FGDB) or shapefile.


texts:
I am trying to analyse the images that I downloaded from MODIS, which are the Vegetation Indices 16-Day L3 Global 1km. With this data I want to analyse the change in vegetation over the years in a area which now suffers with drought.

So, what is the best tool to do this? The files are all in .hdf format, and I do not know how to work with this or convert to a more useful format.


texts: How do I import .lyr files to QGIS? ukhab files. want to add them to a project

texts:
I am trying to add these layer files to QGIS but when I go to import they do not appear in the folder as they are saved as .lyr files but this is not an option.
I am sure there is a really basic thing I am doing wrong!
the files are the UKHabs files.


texts: Configure specific GeoServer workspace in NGINX reverse proxy

texts:
I am using GeoServer 2.12 on CentOS. I have almost 15 workspaces in my GeoServer.

Objective : I want to share my specific workspace to my external users (outside users) so that they can use GeoServer WMS service to QGIS.

Steps followed :

We have installed NGINX on our side and configuring with GeoServer.
I have searched on web but not able to find an example, on GeoServer side they mentioned that one way is to 'isolated workspace'. Main part is how do I configure this with NGINX side ?


texts: Use QGIS field calculator raster_value() function inside Graphical Modeler

texts:
I need to use the raster_value() function (from field calculator) inside a graphical modeler workflow.
The function works if the raster is already loaded on QGIS, and then calling the raster name on the function.
But I need to have an input raster on the model (which might not have been preloaded on qgis) and then use the raster_value() over this raster.

Is there a way to map an inputted raster on modeler and then use the raster_value() function?

I'm using Win10 and QGIS 3.10


texts: Three Dimensional Routing Input Location Errors and Questions, ArcGIS10

texts:
I am having a bit of a problem with my three dimensional route solving. I am creating an interior building network, and ideally would be routing from room to room. I created a line feature class, with z values, implemented it into a network dataset. I built a model that was successfully routing in ArcScene, but something has broken. Whenever I click to interactively add a set of stops to find a route the first one returns as successful but the second one always is returned as an error. The routing does not function if the two stops are not placed on the same line segment, or across a break (simple vertex, or another adjoining segment).

I was also wondering if there was a way to potentially add stops by defining it based on an attribute from a point class that is contained in the network, rather than clicking interactively? I do not wish to add an entire point feature class, just specific points from one.

If anyone has any ideas, they would be much appreciated. Thanks.

I am using ArcGIS 10.


texts: Why does QGIS not show my new geometry columns?

texts:
I am trying to connect a few PostGIS Tables to Quantum GIS.

The problem is that the geometry columns QGIS is showing are not the ones from the PostGIS database. QGIS shows me some columns I used to have before. Newer ones I created later are not shown. So I can't display my actual data.

It seems like QGIS is keeping some old information about my tables. Is there any way to refresh PostGIS connections?

Even deleting the old connections and building new ones does not work.

I checked all the tables on the database level. Everything is alright there. So it seems like a QGIS bug...


texts: GeoServer transform POST request to GET?

texts:
I have the following geoprocessing POST request:

<wps:Execute xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" version="1.0.0" service="WPS" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0    http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
<ows:Identifier>gs:BufferFeatureCollection</ows:Identifier>
 <wps:DataInputs>
   <wps:Input>
    <ows:Identifier>features</ows:Identifier>
    <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wfs" method="POST">
    <wps:Body>
      <wfs:GetFeature xmlns:medford="http://medford.opengeo.org" service="WFS" version="1.0.0" outputFormat="GML2">
      <wfs:Query typeName="medford:firestations"/>
      </wfs:GetFeature>
    </wps:Body>
   </wps:Reference>
  </wps:Input>
  <wps:Input>
    <ows:Identifier>distance</ows:Identifier>
    <wps:Data>
     <wps:LiteralData>0.024</wps:LiteralData>
    </wps:Data>
  </wps:Input>
</wps:DataInputs>
<wps:ResponseForm>
  <wps:RawDataOutput mimeType="text/xml; subtype=wfs-collection/1.0">
  <ows:Identifier>result</ows:Identifier>
 </wps:RawDataOutput>
</wps:ResponseForm>
</wps:Execute>

How can I transform it to GET request?


texts: Select by location of doughnut hole polygon

texts:
I have vectorized the Corine Land class raster and want to select certain classes that are "inside" other classes.

How can I select the polygons circled in blue?

All the select by location options do no work since there are no overlaps.

Explanations per the comments:
What you see is the area around Barcelona, Spain. the red outlines urban polygons, consisting of several types of CLC classes that have been polygonized from the original Corine raster. The green polygons are "green" areas in and around the city.
CLC, as a raster dataset, has only one class for every 100*100 meter pixel through the whole of Europe.
What I'm trying to do is select "green" polygons that are completely surrounded by "urban" polygons, and then export them as a separate layer.


texts:
I have data with emoji characters in a field that in ArcMap are shown as emojis in the Attribute Table.

I've tried all the usual encoding suspects in QGIS but haven't managed to get them display correctly.

Does anyone know if and how it can be achieved in QGIS?


texts:
Im using to gdal_translate to create tiled DEM, tif(float32) -> PNG (UInt16).
All the images trun out black.
When I change UInt to byte - everything work as expected.
What am I doing wrong?

import os, gdal

in_path = 'C:\\xxx\\'
input_filename = 'input.tif'

out_path = 'C:\\xxx\\'
output_filename = 'tile_x'

tile_size_x = 8129
tile_size_y = 8129
ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize

xCnt, yCnt = -1, -1

for i in range(0, xsize, tile_size_x):
    xCnt = xCnt+1
    yCnt = -1
    for j in range(0, ysize, tile_size_y):
        yCnt = yCnt+1
        com_string = "gdal_translate -of PNG -ot UInt16 -srcwin " + str(i)+ ", " + str(j) + ", " + str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(xCnt) + "_y" + str(yCnt) + ".png"
        os.system(com_string)


texts: QGIS Georeferencer--result raster far away and distorted

texts:
I've been trying to georeference a raster.The ground control points I choose using 'from map canvas' are quite close to the original points, but the results are far far away and distorted. Could anyone help with it? The screenshots are below:
1. After choosing GCPs and setting up transformation parameters, you can see the red lines indicating the to-be-moved-to positions:

the GCPs in QGIS:

The result raster is distorted and far far away from original location:


texts: WMS layer is displayed with ruined geometry in OpenLayers

texts:
I want to display a shapefile in OpenLayers using MapServer 7.6.
It is shown properly in the browser using the mapfile's link.

But when I try to display it in OpenLayers, the geometry is completely ruined.

I get no error messages.

What is the cause of this and how could I solve this?

Here is my code:

const map=new ol.Map({
    view: new ol.View({
        center: [2294319.8093533223, 5942407.910713993],
        zoom: 1
    }),
    target: "js-map"
})

    var wmsSource = new ol.source.ImageWMS({
        url: 'http://localhost/cgi-bin/mapserv.exe?MAP=/ms4w/apps/GEO/test.map',
        params: {'LAYERS': 'world'},
        serverType: 'mapserver',
        crossOrigin: 'anonymous'
    });

    var wmsLayer = new ol.layer.Image({
        source: wmsSource,
        title: "world"
    });

const overlayGroup=new ol.layer.Group({
    title: "Overlays",
    fold: 'open',
    layers: [wmsLayer]
})
map.addLayer(overlayGroup);

var layerSwitcher=new ol.control.LayerSwitcher({
    groupSelectStyle: 'group'
});
map.addControl(layerSwitcher);

And the mapfile:

MAP
  EXTENT -180 -90 180 90
  IMAGETYPE "png"
  NAME "world"
  SIZE 1381 1025
  STATUS ON
  UNITS METERS

  OUTPUTFORMAT
    NAME "png"
    MIMETYPE "image/png"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    IMAGEMODE RGB
    TRANSPARENT FALSE
  END # OUTPUTFORMAT

  SYMBOL
    NAME "symbol1"
    TYPE ELLIPSE
    FILLED TRUE
    POINTS
      1 1
    END
  END

  PROJECTION
    "init=epsg:4326"
  END # PROJECTION
  LEGEND
    KEYSIZE 20 10
    KEYSPACING 5 5
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    STATUS OFF
  END # LEGEND

  QUERYMAP
    SIZE -1 -1
    STATUS OFF
    STYLE HILITE
  END # QUERYMAP

  SCALEBAR
    INTERVALS 4
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    SIZE 200 3
    STATUS OFF
    UNITS MILES
  END # SCALEBAR

  WEB
    METADATA
      "wms_srs" "EPSG:4326 EPSG:3857"
      "wms_onlineresource"  "http://localhost/cgi-bin/mapserv.exe?MAP=/ms4w/apps/GEO/test.map"
      "wms_enable_request"  "*"
      "mapmanager_version"  "1.1.0.0"
      "ows_enable_request"  "*"
      "wms_feature_info_mime_type"  "text/html"
      "wms_encoding"    "UTF-8"
      "wms_title"   "world"
    END # METADATA
  END # WEB

  LAYER
    DATA "E:/VECTOR/world_countries.shp"
    METADATA
      "wms_include_items"   "all"
      "wms_title"   "world"
    END # METADATA
    NAME "world"
    PROJECTION
      "init=EPSG:4326"
    END # PROJECTION
    STATUS DEFAULT
    TILEITEM "location"
    TYPE POLYGON
    UNITS METERS
    CLASS
      NAME "world"
      STYLE
        COLOR 58 132 59
        OUTLINECOLOR 30 30 30
      END # STYLE
    END # CLASS
  END # LAYER

END # MAP


texts: Using values from one field as parameters and another field to label in ArcMap?

texts:
I have this attribute table from a shapefile of a road network. I want to labels to display "ROAD_NAME" but I want them to appear at different scales. I would like to use "Length_Cat" to group which roads appear at which scale.

is there a way I can use an SQL expression?


texts: Reducing a feature collection over features instead of columns

texts:
I created a feature collection where features are grid cells, columns are monthly images and the values represent pixel area. I would like to find the mean for each feature over all columns. I tried to use reduceColumns but I get the mean for each column rather than for each feature.

script:

//----load image collection, filter, and create an image with multiple bands

var datasetAll = ee.ImageCollection(MODIS/006/MCD64A1)
                  .filter(ee.Filter.date('2000-01-01', '2019-12-31'))
                  .filterBounds(border);

var DSall1B = datasetAll.select('BurnDate')

var DSallImage = DSall1B.toBands()

 //----filtering grid to the area

 var grid1degreeF = grid1degree.filterBounds(border);

//----creating a layer where pixels value are the pixel area in m^2

var DSallValue1 = DSallImage.where(DSallImage, 1)

var DSallAPP1 = DSallValue1.multiply(ee.Image.pixelArea())

//----summing area in each grid cell over all bands

var DSallAPP1d1 = DSallAPP1.reduceRegions({
  reducer: ee.Reducer.sum(),
  collection: grid1degreeF,
  scale: 500,
  });

Here are my assets:

border- https://code.earthengine.google.com/?asset=users/stavo/amazon_border_harvard

grid- https://code.earthengine.google.com/?asset=users/stavo/Vectorgrid_1degree

This is where I got so far. I end up with the feature collection I described before and now I wish to get the mean for each cell over all columns (dates).

This is the method I tried:

var properties = DSallAPP1.bandNames();

var means = DSallAPP1d1
.filter(ee.Filter.notNull(properties))
.reduceColumns({
  reducer: ee.Reducer.mean().repeat(229),
  selectors: properties
});

After calculating the means I was planning to use reduceToImage so I can end up with an image of the area means over time divided into grid cells. I assume that my method is not the most efficient one. Please advise how to calculate the means or how to approach the issue differently.


texts: Exporting Night Lights data as per the given shapefile in Google Earth Engine

texts:
I am new to GIS and Google Earth Engine Platform. I am trying to get the night lights for all counties in a country. I have uploaded the shapefile as an asset to the platform. I want to export the avg_vis for the region as per the boundaries defined by the shapefile for each county. I want to export them hopefully as an excel file. Additionally, I want to do this for 35 years for which data is available. The code I have till now is:

var shp = ee.FeatureCollection(table);
var dataset=imageCollection.filter(ee.Filter.date('1992-01-01', '1992-12-31')).map(function(image){return image.clip(shp)});
var nighttimeLights = dataset.select('avg_vis');
var nighttimeLightsVis = {
  min: 3.0,
  max: 63.0,
};
var nl=dataset.reduceRegions({

collection: shp,
  reducer:ee.Reducer.mean()

});

    print(table.size());
    Map.addLayer(nighttimeLights, nighttimeLightsVis, 'Nighttime Lights');

imageCollection is the DMSP Nightlights Image Collection and table is the shapefile that I have. I am getting an error that reduceRegions function does not exist. I understand that it maybe happening due to casting but how do I resolve it and also, how to export it subsequently?


texts: 3D roof polygon creation from polylines in ArcGIS for Desktop?

texts:
I am Working in ArcGIS 10.1 for Desktop, and I would like to create roof polygons from lines, the lines are captured with photogrammetry technique.

As per my requirement I have found tool "3D Building Generator Tool", but tool is based on VBA language.

Can anyone guide me in arcgis 10.1 to create Buildings roof polygon using with Arctoolbox or is there any alternative process?


texts: Error while importing shapefile to PostGIS

texts:
I'm trying to import ArcGIS shapefile to PostGIS. I get an error like this:

My client and server encoding both of them "UTF-8"

I get an error especially in the point type.

How can I fix it?

I tried Latin1 and Latin3 but did not work.


texts: QGIS 3 and Python import qgis.core error

texts:
I am trying to set up QGIS3 and python so I can write a script that uses GDAL- I want to do this from the OSGEO4W shell. I have read a lot of questions and followed some tutorials which led me to do this once installing the software:

Firstly I load this cmd which works without error in the OSGEO4W shell:

@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat

@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\Qt5\bin
path %PATH%;%OSGEO4W_ROOT%\apps\Python36\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python36
set PATH=C:\Users\Issy\Documents\QGIS;%PATH%

Then there is the issue:

python3

Python 3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12) [MSC v.1900 64 bit
(AMD64)] on win32

import qgis.core

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\OSGeo4W64\apps\qgis\python\qgis\core\__init__.py", line 27, in ><module>
    from qgis._core import *
ImportError: DLL load failed: The specified procedure could not be found.

I have tried numerous pythonpath settings and installed/re-installed with different settings but I just can't get this to work.


texts: Setting up Raster Calculater for QGIS Python Console: ImportError?

texts:
I'm trying to follow the instructions here How to evaluate raster calculator expressions from the console? to set up the python console to to use the raster calculator.

I tried:

from qgis.analysis import  QgsRasterCalculator, QgsRasterCalculatorEntry

But that gave

Traceback (most recent call last):
  File "<input>", line 1, in <module>
ImportError: cannot import name QgsRasterCalculator

The following seems to work fine, but doesn't get me the raster calculator commands.

from qgis.analysis import *

What is going on, and how do I fix it?


texts: Are LiDAR datasets with same vertical datum, but based on different geoids comparable?

texts:
I have two sets of LiDAR for the same region that I would like to compare. Both are in NAVD88. The older LiDAR set used GEOID03 for the transformation. The newer LiDAR set used GEOID09 for the transformation.

To accurately compare these two LiDAR datasets do I need to update the old dataset to use the new geoid?

The following ngs.noaa.gov link is relevant, but it is still a little unclear to me because in the first sentence it says use GEOID03 if it is from 2003, and then, goes on to say GEOID09 will be more accurate for data use. The gist sounds like that it should be converted to the same geoid. That is in line with my original thinking that it is a systematic error that should be corrected. I was hoping there would be someone out there that had experience with geoid transformations to confirm as it seems like it might add additional uncertainty.

NOAA's quote:

  GEOID03 was fit to the data available in 2003. If you are still using data from 2003, the GEOID03 will fit better. However, NGS has changed the heights in the database for many of these points. Hence, GEOID03 will no longer fit the current values to (in some cases) better than the dm-level. GEOID09 was developed using current values. Hence, GEOID09 should be used now if you want to get results consistent with those showing on the NGS datasheets (our current best estimate of actual heights).


texts: Attaching ArcGIS image server dataset in QGIS?

texts:
I have tried adding the dataset in QGIS using
'qgis.utils.iface.addRasterLayer()'. But it throws this error.

  Invalid Layer: GDAL provider Cannot open GDAL dataset
  http://precisionedgeag.com:8080/arcgis/rest/services/ImageServices/RapidEye_R_NDVI?f=json&pretty=true:
  `C:\DOCUME~1\15937.DSM\LOCALS~1\Temp\file.dat' not recognized as a
  supported file format. Raster layer Provider is not valid (provider:
  gdal, URI:
  http://precisionedgeag.com:8080/arcgis/rest/services/ImageServices/RapidEye_R_NDVI?f=json&pretty=true

How to add this ArcGIS server (Imageservice dataset) in QGIS?


texts: How to tame massive tiled rasters in leaflet webmap

texts:
I have a project that has grown beyond it's initial scope.   I have multiple large rasters tiled in a leaflet webmap.   They were tiled using gdal2tiles.py which produced ~500,000 tiles per raster.

Now I would like to migrate to another server and it's completely unwieldy.   Could a geoserver handle these massive rasters?   What other options do I have to serve large rasters in a leaflet webmap?


texts: Temporal spectral analysis in order to differentiate agricultural/croplands from natural/semi-natural grasslands using remotely sensed data

texts:
I am working on a critically endangered bustard species in India (Houbaropsis bengalensis). I need to prepare a land cover map of its 'suitable' grassland habitat across its distributional range. I have points (and can construct polygons from them using Google Earth's help) where 'suitable' (for the species) grasslands exist since I did presence/not detected surveys in its distribution range.

Since the very same region also has a great amount of cultivation activity, sometimes land cover maps resulting just from the spectral signature using supervised classification tends to pick up croplands as 'grasslands'.

In order to avoid this obviously erroneous assessment, we (myself and collaborators) thought that using a temporal-spectral signature might be a better thing to do since croplands experience this sudden dip in NDVI when harvested and till new crops are sown (they may also have multiple cropping within a year), whereas grasslands do not (well...grasslands too are burnt once a year as a habitat management practice).

May I request people on this forum to help me in making a more detailed land cover map of this study area(s)/distribution range using temporal-spectral signatures in a supervised/unsupervised classification framework in order to distinguish between otherwise similar land cover types such as that in my case (natural/semi-natural grassland and cropland)? Is it even possible to do this in QGIS entirely, or perhaps using the interface of R and QGIS?

Please note that I would be using PROBA-V NDVI (333metre resolution) data for this.


texts: Location estimation methods based on multiple surfaces

texts:
I’m probably reinventing the wheel here and wonder if somebody can point me in a better direction.

I have about 10 surfaces of signal strength (z) received at a transmitter. Here I’ll provide an example with four surfaces r1, r2, r3, r4.

library(raster)
library(gstat)
set.seed(3178)
n <- 100
r0 <- raster(nrows=n, ncols=n,xmn=1, xmx=n,ymn=1,ymx=n)
r0 <- as(r0,"SpatialGrid")

pts <- data.frame(x=c(10,20,60,80,1,1,n,n),y=c(10,80,70,40,1,n,1,n),
                  id=c(1:4,rep(0,4)), z=c(c(rep(1,4)),rep(0,4)))

pts <- SpatialPointsDataFrame(pts[,1:2],data=pts[,3:4])

r1 <- idw(z~1,pts[c(1,5:8),],newdata=r0)
r2 <- idw(z~1,pts[c(2,5:8),],newdata=r0)
r3 <- idw(z~1,pts[c(3,5:8),],newdata=r0)
r4 <- idw(z~1,pts[c(4,5:8),],newdata=r0)

r1 <- raster(r1)
r2 <- raster(r2)
r3 <- raster(r3)
r4 <- raster(r4)

I also have a series of points with signal strengths that correspond to each surface. The locations of these points is unknown. The task is to determine a location for each point. E.g., for unknownXY:

# Extract a point. For the real data I would have the z data but not x and y
cell2get <- floor(runif(1,1,n^2))
unknownXY <- c(extract(brick(r1,r2,r3,r4),cell2get))

In this case there is a four-way intersection of the contour lines that I could find suing some kind of intersection technique (e.g., rgeos::gIntersection).

# Here is the the contour for the unknown point on each raster
plot(r1,axes=F)
contour(r1,levels=unknownXY[1],add=TRUE)

plot(r2,axes=F)
contour(r2,levels=unknownXY[2],add=TRUE)

plot(r3,axes=F)
contour(r3,levels=unknownXY[3],add=TRUE)

plot(r4,axes=F)
contour(r4,levels=unknownXY[4],add=TRUE)

# here are all 4 contours
plot(r0,axes=F,col="white")
contour(r1,levels=unknownXY[1],add=TRUE,drawlabels=FALSE,col="red")
contour(r2,levels=unknownXY[2],add=TRUE,drawlabels=FALSE,col="blue")
contour(r3,levels=unknownXY[3],add=TRUE,drawlabels=FALSE,col="purple")
contour(r4,levels=unknownXY[4],add=TRUE,drawlabels=FALSE,col="green")

But there are times when, I think, that there will not be a completely pure intersection (e.g., the contours will not exactly overlap). In that case I’m wondering about using a probabilistic method for finding x,y for a given vector of z.

One thing I was thinking was calculating distances. E.g.,

# calc distances
r1d <- sqrt((unknownXY[1] - r1)^2)
r2d <- sqrt((unknownXY[2] - r2)^2)
r3d <- sqrt((unknownXY[3] - r3)^2)
r4d <- sqrt((unknownXY[4] - r4)^2)
plot(brick(r1d,r2d,r3d,r4d))

And them summing them:

# Straight sum? Should they be weighted in some manner?
rdSum <- sum(brick(r1d,r2d,r3d,r4d))
plot(rdSum)

I could find the min value to get the ``best'' point:

# Use min to get possible location. Seems easy to get trapped.
rdSumMin <- rdSum == minValue(rdSum)
plot(rdSumMin)

Or use a threshold value to get a series of locations that are the most likely:

# Or use quantile to get possible locations?
thresh <- quantile(rdSum, probs = c(0.01))
rdSum99 <- rdSum < thresh
plot(rdSum99)

Is there a way of finding a surface that is explicitly probabilistic?

What I’d like to do is be able to determine possible location with x% likelihood.


texts: Which algorithm is used for fill sinks in arc hydro tools?

texts: How I can get information regarding the % of voids filled using fill sinks in DEM and the algorithm used to fill the voids using arc hydro tools?


texts: I am working on a manuscript using the gridded soil type data product obtained from worldgrids.org and need to cite the data source. However the domain has expired. Does anyone have any idea what happened, or if the datasets formerly on worldgrids.org are hosted elsewhere?


texts: Adding polygon to layer from text file using QGIS

texts:
I am completely new to QGIS and to GIS in general. I have start to work my way through the documentation but have a question that I cannot find a answer to. When looking into the documentation related to Data Capturing I find plenty of examples when creating lines and polygons using interactive picking.
If I instead have the data in text file (containing coordinates) how do I import the text file?
I have tried to create a layer from "delimited text files" but I do not want to create a new layer I would like to add this data to an existing layer and in addition add attribute data.


texts: How to get a more refined and smooth DEM?

texts:
I have Two DEMs (Digital Elevation models)

DEM (Image 1) - received from my institute
DEM (Image 2) - downloaded from srtm

Both DEMs have been styled similarly. But DEM (Image 1) looks better than DEM (Image 2). The resolution of both the DEMs is also same that is 30m X 30m.
My question is why they look differently after doing same styling in qgis and how I can get such DEMs as in image 1.
My OneDrive link to download the images is
https://1drv.ms/u/s!At3AD05XRRsNiFQRmp2KCiO1Mkmf?e=i6Owf9


texts: Missing features after joining by attribute, but geometry is valid

texts:
I've tried joining two point layers by location on QGIS 2.18.27 (Mac iOS 10.13.5 High Sierra), but some points in the new vector don't end up appearing. One point layer is an interpolated point layer of a line, with one point every 500 meters (we'll call this interpolated_point layer), and the other point layer is a point at every node of above-mentioned line layer (we'll call this the node_point layer), both set to World Mercator [EPSG:54004].

Points and attribute table for the node_point layer

Points and attribute table for the interpolated_point layer

The node_point layer (left) and the new joined layer (right, using the Join attributes by location tool) with missing feature highlighted

I used the Join attributes by location tool in the Vector > Data management tools menu and set the node_point as the Target layer and the interpolated_point layer as the Join vector layer. The Geometric predicate is within and Precision is set to 500 m (so that if a point on the interpolated_point layer isn't within another point in the node_point layer, it'll try to get the attributes of the point within 500 meters; that is, that's how I think it works).

Join attributes by location dialogue box with entered parameters. points_mercator is node_point and interPoints_mercator is interpolated_point

I've tried checking the validity of the node_point layer and all the features node_point layer are all right, so I'm not sure what I'm doing. Reversing the Target and Join vector layer parameters don't change the (undesired) result either. Any thoughts? The expected result is that I have a new layer that has all and only the points from the node_point layer with the attribute tables combined by location.

Again, using QGIS 2.18.27 on a Mac iOS 10.13.5 High Sierra.

EDITS:

Added screenshot of interpolated_point and node_point layers and their attribute tables (as suggested by @pdavis)
Added screenshot of Join attributes by locationdialogue box with parameters used (as suggested by @pdavis)
Added expected results (as suggested by @pdavis)
Added CRS used


texts: QGIS 3 How to Create Glow on Layer with Gradient Symbology

texts:
I have a line layer (street network) coloured with a gradient scale based on certain attributes. I would like to add a glow effect from the draw effects properties option, however the glow can only be a solid colour or a gradient in the glow.

What I would like is that the glow matches the colour of the line segment. For example, some lines are green and I would like them to have a green glow. Other lines are blue and I want them to have a blue glow. I can't figure this out! Any suggestions?


texts: Openlayers 3 - how to get id/index ol.geom.LineString vertex by coordinates

texts:
I'm having a list of points (coordinates + extra info). I want to show them on the map as LineString and when I hover a vertex/point I want to get a callback where I display extra info for the given point. As the list could be big I don't want to create a feature for each point, but just one for the LineString. Otherwise I end up in performance problems.
I can easily take the closest point from the LineString, but it's just coordinates, without any identification ... index would be sufficient.

Here is how my code looks like:

var geometry = new ol.geom.LineString(positions);
map.on("pointermove", function(event)
{
    if (event.dragging) { return; }

    var mouseCoordinates = map.getEventCoordinate(event.originalEvent);
    var closestPoint = geometry.getClosestPoint(mouseCoordinates);
    console.log(closestPoint);
});


texts:
I have a WMS layer hosted on Geoserver that is being shown on a web map using the Google Maps API.

The problem is that for now I want to use basic authentication on the Geoserver. Every time you try to access the server (or the layer) in browser, it will ask for a username and password. (Note: The prompt only shows up on Firefox. Chrome just logs 401 errors in the console).

What I've been trying to do is authorize it in the background so the user does not have to enter any credentials. Here's what I tried.

Before loading the map, send a request with credentials to fetch at least one tile: http://user:password@localhost:8080/geoserver/wms?layers=mosaic&version=1.1.1&Service=WMS&request=GetMap&Styles=&format=image%2Fpng&SRS=EPSG:4326&BBOX[my bounds]=&width=256&height=256. My server-side code is in Python so I used the requests library on that URL to at least get something about the layer, like status code or content type. This works fine.
My assumption was that this would have authenticated the access to this layer (and server) for the user in the background.
But when I load my WMS layer on my map, it will ask for authentication again on Firefox, and Chrome just logs 401 errors. The only difference is that in my JavaScript file, the URL is without user:password. I don't want to expose the credentials on the client-side.

What else can I try to work around this issue?


texts: LandXML to TIN places TIN in wrong location

texts:
I'm using ArcMap 10.8 with 3D Analyst. I received a LandXML from a civil engineer, who I assume used Civil3D or similar to create it. The LandXML has several TINs representing existing surface and proposed surfaces of a proposed project. When I use the LandXML to TIN tool, the TIN is placed in the incorrect spatial location.
You can see in this screenshot the overpass is shifted to the northeast by about 575 feet.

I have tried the following:

Apparently the data uses a modified state plane with a scale factor. I have tried editing the xml with the scale factor before converting to TIN, but the result is the same.
I have tried editing the prj.adf file in the TIN folder with the scale factor applied, same result.
The TIN result also has no projection information, despite specifying in the Environment settings. I set the XY coordinate system in the properties of the TIN to the state plane as well as modified state plane with scale factor applied, both of which have no effect on the location of the TIN.
I requested a dwg from the same engineer so I could inspect it; this data lands in the correct location when added to ArcMap.
Using TIN to Raster, I converted one of the proposed surface TINs as well as the existing surface TIN to .tif. I georeferenced the existing surface raster to an approximate location based on aerial imagery (getting desperate). I copied that .tfw file and put it in the same folder with the proposed surface raster, and renamed it with the name of the proposed surface raster. This did not work the same way as a .wld file would with CAD, and resulted in a shrinking and misalignment of the proposed surface raster.

How do I get this TIN to land in the correct location?


texts:
I have problem with QGIS (2.0.1), PostGIS (2.0.4) combination.
- A can load tables or views by Add PostGIS layer. This is OK.
- But a can't load data by DB Manager and SQL Window. Executing query run but load result is problem.

Error message is:

  PostgreSQL layer has no primary key

But Primary key exist in table. DB Manager find out them itself.

Do you have any idea where is problem?
Could my problem be with different schema then public?

Edit: Do this when I click to Load Now only. And realy do this outside public schema. In public schema it's OK.


texts: multiple glitches in rendered rasterfile after reprojecting in QGIS

texts:
I reprojected a natural earth raster file (Natural Earth II with Shaded Relief, Water, and Drainages, CRS WGS84) with a costum CRS (azimuthal orthografic) using Raster->Projections->Warp in QGIS. Worked fine except glitches at the international date line and at the horizon of the hemisphere (see attached image). An idea how I can prevent these?

Edit: as @user30184 suggests adding -wo SAMPLE_STEPS=100 and -wo SOURCE_EXTRA=10 fixes the glitches


texts: Is there an R function for creating spatial weight matrix for “spatial difference in differences” method?

texts:
I am trying to apply "spatial difference in differences" method in R for a US county's block groups (n block groups) over 2 different times (2000 and 2015).
In order to do this spatial regression over two years, I need to create a spatial weight matrix based on a contiguity list of neighbors.
Based on a paper by Delgado(2015), this spatial weight matrix is a

(2n*2n) block-diagonal row-standardized spatial weights matrix containing non-zero elements for spatial units belonging to contiguity class.

Can anyone let me know of an R function for creating such a neighbors' weight matrix that can be created for a county block groups' shapefiles over two years and then used for spatial difference in differences specification?


texts: Creating Polygons Based on Point Count in ArcGIS for Desktop?

texts:
I am looking to create polygons based on a point count, and hopefully find a way to automate this too. My example is, I have 1,000 points in my city boundary and I want to create polygons over the area of 100 people in the city boundary. So my result would be 10 polygons each covering a 100 people and not overlapping. I am using ArcGIS Desktop 10.3 with an Advanced level license.

I couldn't find the answer online.

I would want to automate this too so I could do this multiple times as the population grows.

I don't want to go by distance I want to go by point count alone.


texts: Seeking polygon shapefile of countries, states and islands?

texts:
As often happens, the hardest part of GIS is to come up with clean shapefile with an attribute table.

I need to create a choropleth map, but all I have is this:

Domicile,Captives
New Zealand,24
Vanuatu,7
Hong Kong,3
Labuan,38
Micronesia,14
Singapore,64
British Columbia,21
Anguilla,252

As you can see, some of them are countries, some are US states and there are islands too.

The way I think of doing this is to first look for a world polygon shapefile of countries, then a shapefile of states and then look manually for what is left.

Is  there a shapefile that contains everything?


texts: graphics all change to same color when map zooms in or out

texts: As the title states, I have a working "draw tool" where you can draw on the map using polygons, simple marker symbols, lines, and text. Whenever you draw any of these symbols on the map using different colors and you zoom in or out all of the symbols will change to the same color.


texts:
I'm trying to strip the last character in a string if when the string is split by spaces the last element has a length of 1. I then want to put that letter in front of the rest of the other elements.

So, if I had a value of 1st Street W, I would like to convert it to W 1st Street

I'm using this code in the "Pre-Logic Script Code" box in ArcMap field calculator:

def getPrefix(Street_Name):
  if len(!Street_Name!.split(" ")[-1])==1:
    !StreetName2! = !Street_Name!.split(" ")[-1] + " " + !Street_Name!.split(" ")[0:-1]

And in the box below that:

def getPrefix(!StreetName2!)

I'm getting an error 000989, syntax indention error, but I have everything at two spaces in cascading from calling the function. Pretty sure that's the ESRI specification for the field calculator.


texts: How to find nearest neighbors of polylines?

texts:
I have a shapefile containing irregularly placed polylines (i.e., survey routes). Following the general idea of a paper I read, I would like to use tessellation to identify the neighbors for each polyline, but I have no idea how to accomplish that.  I see there is a tool to create Thiessen polygons around points which (I think) is similar to what I want to do, but I don't see a similar option for lines.  I am using ArcGIS 10.1. Any specific ideas on how to proceed would be much appreciated!

Additional information

1) The paper is here:  http://www.esajournals.org/doi/abs/10.1890/03-5247. Their brief description: "We delineated a spatial neighborhood on an irregular lattice by tessellating the sample routes, and from this neighborhood structure derived an adjacency matrix of first-order neighbors."

2) The routes do not overlap.

3) I do have access to Spatial Analyst.

4) The goal is to determine the nearest neighbors of each polyline so we can control for spatial auto-correlation of count data from the routes.


texts:
I am trying to draw a line (vector layer) of a certain length that that tries as much as possible to dodge the green polygons.

Assuming a starting point, is there some sort of line centroid algorithm/function already created that I could use?


texts: Why is "Join Attributes by Location" reporting different CRS (and running forever)

texts:
I have a point feature layer with ~240,000 parcel centroids. I have a second feature layer containing ~9000 subdivision polygons. I want to assign the name of the associated subdivision as a new field to each point in the parcel centroid layer. The polygon layer is a shapefile containing a single field, the subdivision name. The point layer is a csv file containing 3 fields - the parcel ID, lat, lon and is loaded using layer->Add Delimited Text Layer.

The polygon layer is in WGS84/EPSG:4326. I have explicitly set the point layer to the same CRS. Everything looks sane on the map. I then try to do a spatial join: Vector->Data Management Tools->Join Attributes by Location. I set the target layer to the point layer, the join vector layer to the polygon layer. I set attribute summary to "Take attributes of first feature". I set the output shape file and click OK. I immediately get a warning message that the two layers have different CRS. But I have verified that are both EPSG:4326, as is the project CRS. I click OK and the progress bar immediately shows 15% and then life comes to a halt. The task monitor shows qgis running at 100% of CPU, but after several hours, it's still running at 100%, and the progress bar has not updated. I spatially indexed both layers.

Why did I get the different CRS message?

Why is this taking so long to run? (qgis 2.0.1 Dufour, Windows 7/64bit, Intel i7 2.4GHz)

Is there a better/different approach I should try?

Edit:
For what it's worth, ArcMap 10.2 did the spatial join in about 15 minutes (rough estimate - I wasn't actually timing it.)

If anyone wants to try this join, the data can be found here


texts: Create multiple "Values" from ModelBuilder Iterator

texts:
I am using "Iterate Row Selection" in ModelBuilder to run through every line of a table. I am able to have the iterator create a rolling variable (Value), by selecting a field in "Group by Fields."

If I "Group by Fields" on the field CNTRY_NAME, the iterator creates a value for that record on each iteration (i.e. Afghanistan). I can then plug the variable into a different part of the model as an input.

I would like it to create a second value for an additional field (i.e.  CNTRY_ID) on each iteration, based on the current selected record. By default, the iterator only creates one value.

If I add an additional field into the Group By Fields, the fields are concatenated, and I still only get one "Value" variable.

Is there a way to create an additional "Value" variable from an iterator in ModelBuilder?


texts: Losing DEM (from LAS Dataset) data coverage when calculated in higher resolution?

texts:
I am trying to calculate DEMs for different resolutions (15 x 15, 7.5 x 7.5 and 3.75 x 3.75 m) using LAS Dataset To Raster tool in ArcGIS 10.5.

Unfortunately, no matter which settings I use for interpolation, the output DEMs do not cover the whole Las Dataset area (see image):

It seems that the calculated DEMs have shrunk proportionally the same in both (x and y) directions. The LAS Dataset extent is 27 x 9 km; however, if I use another LAS Dataset (with larger or smaller extent) I have the same issue.

I assume it is something in the cashing memory of the tool, or the environment settings, although I haven't changed anything there so they should be default.

Does anyone know what is going on there?


texts: Mapbox Studio best practices for line text-placement? Can you force 1 label per segment?

texts:
I'm struggling with text-placement: line; in Mapbox Studio. I'm fighting between 1 label being cut-off and 2 labels cutting each other off.

I'm trying to label fjords in Alaska. Originally, I had the line that followed the length of the fjords, but when I labeled them in Mapbox Studio the label repeated again and again. The fjords make a large Y shape so I broke up the lines into 3 legs roughly the length of the text I want to appear.

In one section the label goes from being cut-off to 2 labels cutting each other off.

I've tried tweaking 'text-character-spacing' and 'text-size' with no luck and I'm tempted to edit the lines in some way of getting the length just right, but it could take hours to tweak everything just right.

Easiest solution I think would be if you could force 1 label per segment but I can't find any documentation or examples of doing that.

Any suggestions?

Edit: I need to mention I've also tried 'text-repeat-distance' and 'text-min-distance' and they don't seem to do anything.

Edit 2: I created a gist in case someone wanted to work with the data and test it out. Gist with link to data source and cartocss


texts: I am in the process of developing distribution map in MAXENT. My species habitat is subset of another mammal habitat. I have distribution layer digitised for mammal species, rasterised and converted to asc. Apparently the raster dimensions are slightly different. If I edit the asc file in word the features disappeared. Need help on how to go about it.


texts: Factors in computing scoring for place name searches, like Geonames

texts:
As part of the FULL response style, the GeoNames search API returns results like the following:

<geonames style="FULL">
  <totalResultsCount>510</totalResultsCount>
    <geoname>
    <toponymName>Las Vegas</toponymName>
    <name>Vegas</name>
    <lat>36.17497</lat>
    <lng>-115.13722</lng>
    <geonameId>5506956</geonameId>
    <!-- abbreviated -->
    <score>92.10135650634766</score>
  </geoname>
</geonames>

How is the score element calculated? I'm guessing there's some sort of search index (e.g. something running on top of Lucene, like Solr or Elasticsearch).

I'm less concerned about exactly understanding the scoring algorithm they use, but rather an understanding which factors would be the most useful in implementing a scoring mechanism if I were to implement a similar web service locally.


texts:
I am having problems understanding how you can ordering the output of a sapply list based on the column headings.  I want to calculate the % of each raster type within a polygon.  This is almost identical to the example/problem from Curlew in Get Raster Values from a Polygon Overlay in Opensource GIS Solutions and thecrashlandingdodo in Step-by-step: How do I extract Raster values from Polygon overlay with Q-GIS or R? but I am experiencing an issue with a sapply output.  Essentially I can't "fix" the raster value % to the correct column in my summarised table, issue shown below.

I am using R 3.2.3 x 64 via RStudio on a windows 7 machine.

My code is

require("raster")
require("rgdal")

# This is a simple 1m (grid) raster with 4 numerical types 1, 2, 3, 5.
Rast <- raster('SpeciesType.tif')

# polygon data same extent as raster dataset
Vec <- readOGR(dsn="K:/Data", layer="StudyAreas")

extra <- extract(Rast, Vec)

extra.count <- sapply(extra, table)

# The % of each raster numerical value within the polygon area
extra.list <- sapply(extra.count, FUN=function(x){ x / sum(x) } )

extra.percent <- as.data.frame(t(sapply(extra.list,'[',1:length(unique(Rast)))))

So the first list (extra.list) is correct I.e. % of each numerical raster type in a polygon e.g.

head(extra.list)

[[1]]

         1          2          5

0.30587297 0.02657454 0.66755248

[[2]]

          1           5

0.007246377 0.992753623

[[3]]

5

1

[[4]]

5

1

[[5]]

        1         5

0.2806617 0.7193383

[[6]]

5

1

But the sapply list does not put the raster value % into the correct column headings e.g. the % in polygon of numerical value 5 is in column 5, then 2, then 1 in the first three rows.

head(extra.percent)

            1          2         5 NA

1 0.305872974 0.02657454 0.6675525 NA

2 0.007246377 0.99275362        NA NA

3 1.000000000         NA        NA NA

4 1.000000000         NA        NA NA

5 0.280661723 0.71933828        NA NA

6 1.000000000         NA        NA NA

How can I alter the data frame/sapply statement to fix this?

extra.percent <- as.data.frame(t(sapply(extra.list,'[',1:length(unique(Rast)))))


texts:
I have OpenStreetMap (OSM) data that I've loaded into a Spatialite database.

How do I get it to look like OSM?

I shouldn't have to add and style each layer individually but I can't figure out what to do next.


texts: GDAL Java Bindings .SetWellKnownGeogCS error

texts:
I am using GDAL 2.3.0 and Java bindings from http://www.gisinternals.com/query.html?content=filelist&file=release-1911-x64-gdal-2-3-0-mapserver-7-0-7.zip

I would like to create a SpatialReference from EPSG code.
This works fine:

hProj = new SpatialReference();
hProj.ImportFromEPSG(32632);

but

hProj = new SpatialReference();
hProj.SetWellKnownGeogCS("EPSG:32632");

generates

Exception in thread "main" java.lang.RuntimeException: OGR Error: General Error
    at org.gdal.osr.osrJNI.SpatialReference_SetWellKnownGeogCS(Native Method)
    at org.gdal.osr.SpatialReference.SetWellKnownGeogCS(SpatialReference.java:473)
    at rgeo.lesson_x.JavaGDALTests.main(JavaGDALTests.java:22)

What could be the reason and how to fix that?


texts: Field calculator error and Geoprocessing results window empty

texts:
I need to use the field caculator for a new field that I have just added but I keep getting the following error:

The problem is that when I open the Geoprocessing Results window is just completly empty, which is so weird, I don't know what to do and I need to use the field calculator to keep going with my work

The field is a short integer and is just equal to another field called OBJECTID


texts: Definitive list of QGIS multiprocessing tools?

texts:
This question is not as straightforward as it seems, but can be summed up with the following:

Is there a definitive list of geoprocessing tools for QGIS that currently support multithreaded processing? The intent is to use more than a single core.

It is well understood that multithreaded rendering is available in QGIS 3, but this has no effect on geoprocessing operations, and in the question history on the various stacks seems to commonly be confused by respondants.

Background:
Datasets are becoming more detailed as acquisition methods increase in detail and density. The requirement for detailed analysis combining multiple datasets and myriad steps to process means that multithreaded processing is an attractive option. For example, ESRI has introduced a small but important handful of multithreaded processing tools in ArcGIS Pro - Pairwise Intersect, Pairwise Buffer, and Pairwise Dissolve. These tools have introduced significant speed enhancements to the traditional intersect, buffer, and dissolve functions when handling large datasets.

QGIS interfaces with a number of GIS toolchains - GDAL, GEOS, Python3 / Cython / Shapely / PostGIS - each which has a different position on whether they support multithreaded or parallel processing. It's really very difficult to know whether a tool you wish to use in QGIS supports multithreaded processing, or which tool from which toolchain might offer a speed improvement without a list of tools that support multithreading. This becomes even more important to users of large datasets when building models that might need to be run with semi regularity.

Despite extensive searching I'm coming up with very little.

So, to answer this question are you either able to point to a preexisting list, or your own list of multithreaded processing options for QGIS?


texts: How to use GRASS raster tools in the Processing modeler?

texts:
In the Graphical Modeler of QGIS, I don't see a way to import a layer into GRASS. The modules "import loaded raster" are not accessible in the Modeler.

I want to use GRASS reclassify module, because SAGA's doesn't save the data entered as a "simple table", forcing me to re-enter the values in the table each time I want to run the model.

But GRASS reclass modules don't work in the Modeler for a lack of loaded rasters.

The rasters I use in the modeler are first treated using non-GRASS algorithms. So I can't just load the rasters into GRASS outside of the Modeler.

Here I may be wrong. Could I load all rasters into the GRASS plugin before opening the Modeler? And then use the GRASS algorithm via the Modeler on them?


texts:
Why doesn't my SLD validate?

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/sld
http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd" version="1.0.0">
<UserLayer>
  <Name>raster_layer</Name>
  <UserStyle>
    <Name>PM 2.5</Name>
    <Title>PM 2.5 Layers</Title>
    <Abstract>Color Scheme for PM 2.5 Rasters</Abstract>
    <FeatureTypeStyle>
          <FeatureTypeName>Feature</FeatureTypeName>
      <Rule>
        <RasterSymbolizer>
            <ColorMap type="intervals">
              <ColorMapEntry color="#fafafa" opacity="0.0" quantity="0.0" />
              <ColorMapEntry color="#fafafa" opacity="0.5" quantity="1.0" />
              <ColorMapEntry color="#fa0808" opacity="0.5" quantity="100.0" />
         </ColorMap>
        </RasterSymbolizer>
      </Rule>
    </FeatureTypeStyle>
  </UserStyle>
</UserLayer>
</StyledLayerDescriptor>


texts:
Arcgis 10.2

I currently have a raster file that I want to turn into a set of polygons. This will help me to create specific buffers around certain points. The extent of my values are 0, and 1. When I change from raster into polygon, I end up with polygons for both 0 AND 1. I only want polygons for pixels with value 1. If I have polygons for values of 0, I end up capturing a bunch of empty space that it totally not necessary to the project.

Here is an image of what I am working on (white = 1, black = 0):


texts: Default field value check in qfield form

texts:
Using Qfield Project, in the input data form is always avaible for all field the check to set the default value for the next features.
It is possible with some settings define if this check will be available or not for each field?


texts: PostGIS Raster ST_AsRaster: error when creating raster from TIN surface

texts:
ERROR:  rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB
********** Error **********

ERROR: rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB
SQL state: XX000

---------------------------------------------------------------------------------

I am trying to create Raster from TIN surface created using number of survey points available to me and store it in database. These points have Z Dimension.

I have created geometry that represents all points using ST_Union and applied ST_DelaunayTriangles to obtain "Delaunay triangulation" around the given input points.

To obtain TIN geometry, I have used workaround provided in below link
https://smathermather.wordpress.com/2013/12/18/2-5d-tins-in-postgis/

I verified the TIN surface in QGIS at this point and it looks fine.

Now I need to create DEM raster for this surface.
When applying ST_AsRaster,

select ST_AsRaster(ST_DelaunayTriangles(ST_Union(ARRAY(select survey_point from table))),150, 150, '2BUI');

I get above error.


texts: Is starting names with numbers a bad data naming convention?

texts:
My company uses ArcGIS and has a project and data file naming standards in place and (for the most part) followed.  Something that has always bothered me about he naming standards is that it mandates starting all project and data file names with the project number - an eight digit number.  I've always held the belief that naming GIS files starting with numbers is a bad thing, and have had (especially with GRIDS) processes fail because of the file name.

I'm looking to amend the corporate standards to drop the project number requirement, however I cannot find much in the way of documentation on why "numbers as a first character" in filename is a bad thing.

Can anyone point me in the right direction as far as resources to support this argument?


texts: Split features script returning Python ERROR 000840: The value is not a Feature Layer

texts:
I am trying to split a feature layer by rows. I've previously some time ago been able to use the below script, however now when I use it, it returns:

  Python ERROR 000840: The value is not a Feature Layer

import arcpy

outputNum = 5000000
outputFCName = "Flows"

def listSplit(myList, n):
    for i in xrange(0, len(myList), n):
        yield myList[i:i + n]

arcpy.env.workspace = r"C:\...\Database.gdb" #The dots are replaced with a proper path
lyr = arcpy.mapping.Layer("Layername")
arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")
fList = list()

with arcpy.da.SearchCursor(lyr, "OID@") as cursor:
    for row in cursor:
        fList.append(row[0])

listGroup = listSplit(fList, outputNum)

for x in listGroup:
    lyr.setSelectionSet("NEW", x)
    arcpy.CopyFeatures_management(lyr, arcpy.CreateUniqueName(outputFCName))

The reason for splitting is that I have a dataset with 24.000.000 rows, which I would prefer in smaller parts.


texts:
I'm a newly minted GIS tech.

At first, using Maplex, the labels in my point file were displaying in what looked like a clumped, truly random pattern: some areas would have labels stacked on top of each other, while other areas would have no labels at all.

My supervisor thinks this had something to do with how the labels are weighted.

Today, the file is displaying all the labels, but they look like this. Please note I do not have "Never Remove" checked.

How can I label my map legibly?

John W. Fell commented:

  From what I can tell it appears that your point cloud is rather dense.
  You may need to produce maps indexed by grid and at larger scales for
  each individual map. Data driven pages may be a good solution.

That's actually what we're doing. This page is a 1:12,000 scale map of a square mile, and is page 11 of 38 in the atlas created by the data driven pages.


texts: GDAL_Translate output .TIF is VERY small compared to map (in QGIS)

texts:
using QGIS 3.2 ... Here are the steps that I followed:

1) Used NCRA (NCO) to generate average dewpoint temperatures for given months in given years (calculating the meteorological seasonal average)
2) Used gdal_translate (command line) to create .TIF files to use in QGIS.

Filename: dewps_fall_1980_1989.nc
Subdataset: dpt
Target PRS: 102004 (USA LCC)

I used the following command-line:

gdal_translate -a_srs EPSG:102004 -of Gtiff NETCDF:dewps_fall_1980_1989.nc:dpt test_new3.tif

It creates the look of the output I want, but it is astronomically smaller than my map despite my PRS settings being correct.

Any ideas? Heres a look at the properties box for the .TIF. Maybe the problem is with the extents? Not sure.

But when I just simply try to export it as a new TIF, correcting the extents to my map (just to try) it always says that it is an unsupported data source and doesn't even create the exported file.


texts: GDAL: gdalwarp wrong transformation to "South_Pole_Lambert_Azimuthal_Equal_Area" EPSG:102020

texts:
I am trying to reproject a raster that is defined in EPSG:4326 to EPSG:102020.

Command used

gdalwarp -ot Float32 -s_srs EPSG:4326 -t_srs EPSG:102020 -r near -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 etopo1.tif etopo1_warped_cli.tif

but the output has a wrong coordinate reference system, as it is shown as

+proj=laea +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

instead of

+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

On the other hand if I manually specify the parameters for the target srs then transformation is ok and the raster is recognized as EPSG:102020 in Desktop GIS software:

gdalwarp -ot Float32 -s_srs EPSG:4326 -t_srs "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r near -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 etopo1.tif etopo1_warped_cli.tif

NOTE:
the above is on Windows, because on Linux the operation returns

ERROR 6: EPSG PCS/GCS code 102020 not found in EPSG support files.  Is this a valid
EPSG coordinate system?
ERROR 1: Translating source or target SRS failed:
EPSG:102020


texts:
I've searched here as best as I could, and can't find a direct answer on how I could go about caching the map to speed up loading on zoom. It seems most caching options were introduced in 1.0.0 which I unfortunately cannot use (I can't leverage it at all per Windy guidelines) A live demo of my project can be seen at https://codespent.github.io/dcsg-weather-map and the codebase at https://github.com/CodeSpent/dcsg-weather-map.

I'm restricted to Leaflet 0.7.7 as I'm using the map provided by the WindyCom API for weather patterns. Is anyone aware of any references they can point me to that'll help with the zoom loading?

I do believe its the Windy overlay that causes so much of the hangup, but if I can at least better some practices to accommodate that'll suffice for me.

Primary use for this will be a 24/7 screen with reloads on a schedule in our operations center for quick reference.


texts:
I work on a QGIS Plugin in python. This plugin import in QGIS project several print layout (A4/A3/Landscape/Portait format) from a directoty of my laptop to the layout manager. The plugin add also a vector layer in my QGIS project. I would like that when I open a layout from the layout manager, the map of this layout is centered on the extent of the layer in my project. I have tried the code below but it does not zoom in correctly on the layer.
project = QgsProject.instance()

zip_layer = QgsVectorLayer("C:/Users/Thysma/Documents/Stage/Data/zip.shp", "Aire d'étude",'ogr')
QgsProject.instance().addMapLayer(zip_layer)

composition = QgsPrintLayout(project)
document = QDomDocument()
template_file = open('C:/Users/Thysma/Documents/Stage/Composeurs/A3 Paysage.qpt')
template_content = template_file.read()
template_file.close()
document.setContent(template_content)
composition.loadFromTemplate(document, QgsReadWriteContext())
project.layoutManager().addLayout(composition)

layout = project.layoutManager().layoutByName('A3 Paysage')
map = layout.itemById("map")
map.setExtent(zip_layer.extent())
map.attemptMove(QgsLayoutPoint(0.150, 21.610, QgsUnitTypes.LayoutMillimeters))
map.attemptResize(QgsLayoutSize(342.235, 275.240,QgsUnitTypes.LayoutMillimeters))

The zoom made gives me the following rendering (the map is not zoomed in on the total extent of the layer):


texts: Flat-ended buffers on touching polygon features in QGIS?

texts:
Please note that emphasis is on POLYGON here NOT POLYLINES!

I need

Buffer/s on one polygon,
with two touching features
and straight touching parts after difference has been applied

I do not want this result:

but, rather something like this (the solid black lines, I've scribbled over the lines I don't want in blue)

I need the buffers to end at the shortest line between the buffer boundary (of the union of the features), and the exterior ring intersection points of the two features.

This solution does not work for Polygons as far as I attempted: How to achieve “Flat-End” line buffers in QGIS?


texts:
I have a points shapefile in QGIS with 40 groups of 20 points each (total of 800 points). For each point I have stored many information.

I need to check this information and add more moving on the map by car.

How can I share the QGIS coordinates of each point to an Android smartphone to get each point?

It would be great to share it with an online service (e.g. importing in Google Maps account or something similar)

ORIGINAL QUESTION WAS: I need to share a 800-point QGIS shapefile with android smartphone for offline navigation.

Considering the high number of points I would need to set up some filter so that navigation is as easy as possible.

I use QGIS 2.18 and different standard android smartphone


texts: How to update PostGIS column with distance from features in the same table

texts:
I have a table 'mytable' with columns 'Geometry', 'featurecode', 'height' and 'pkey'.  There is a spatial index and im using PG9.4

For most of the table the height column is null.  Only where featurecode=2 is there a height value.

I am trying to find an efficient way to update records where featurecode=10 with the height of the closest feature of type 2.  This is a very large table so im thinking I may need to limit the search to items within distance of X meters

Has anyone tried doing something similar?


texts: How to count particular DN value's number in image in Google Earth Engine

texts:
I want to count the the number of particular DN value. I have tried to use mask and loop to do it, but it's so slow if the image is big. Here is my code (the variable 'image' is a image):

    var y =[]
for (var yr= 1856; yr <= 1890; yr = yr + 1) {
 var mask =imgae.eq(yr)
 var sumDictionary = mask.reduceRegion({
   reducer: ee.Reducer.sum(),
   region:aoi.bounds(),
   scale: 30,
   maxPixels: 1e9
 });

  y[yr-1856]=sumDictionary
}


texts: GeoDjango: GDAL library is not available to transform() geometry

texts:
I am trying to take points and buffer them to create polygons in my django project. I queried for points and tried to take the geometry of the first point and transform it to a meters based system (3395) so I can buffer in meters.

In [32]: from django.contrib.gis.geos import Point

In [33]: p = workingqueryset[0].geo_location

In [34]: p
Out[34]: <Point object at 0x1cfeb4bf090>

In [35]: p.transform(3395)
---------------------------------------------------------------------------
GEOSException                             Traceback (most recent call last)
<ipython-input-35-ee3c41c5bae3> in <module>()
----> 1 p.transform(3395)

C:\Anaconda3\lib\site-packages\django\contrib\gis\geos\geometry.py in transform(self, ct, clone)
522
523         if not gdal.HAS_GDAL:
--> 524             raise GEOSException("GDAL library is not available to transform() geometry.")
525
526         if isinstance(ct, gdal.CoordTransform):

GEOSException: GDAL library is not available to transform() geometry.

I have gdal installed with python bindings and even have the GDAL_DATA variable set up for postgis. I've never ran into a gdal issue before. I found this answer but I definitely have gdal already installed and variables set up.

In [36]: import gdal

In [37]: p.transform(3395)
---------------------------------------------------------------------------
GEOSException                             Traceback (most recent call last)
<ipython-input-37-ee3c41c5bae3> in <module>()
----> 1 p.transform(3395)

C:\Anaconda3\lib\site-packages\django\contrib\gis\geos\geometry.py in transform(self, ct, clone)
522
523         if not gdal.HAS_GDAL:
--> 524             raise GEOSException("GDAL library is not available to transform() geometry.")
525
526         if isinstance(ct, gdal.CoordTransform):

GEOSException: GDAL library is not available to transform() geometry.


texts:
I work with ArcGIS 10.3 (Python 2.7.8) and I need to change the grid interval in 70 MXD's.

I tried to find Python code that can do it. Currently, the interval in each MXD is 200 meters and I need to change it to 750 meters.


texts:
I have an 8-bit greyscale image that I want to convert to a "heat map" based on data values (0-255), using a color ramp.  I'm aware that gdaldem works with color-relief, so I used the following steps from a prior post:

Take a raster of your choice (like SRTM as an example)
Create the following text file named col.txt:

0 black
100 blue
200 yellow
300 orange
400 red
500 white

On command line, run

gdaldem color-relief N51E007.hgt col.txt out.tif

When I do this, I get the following:

"ERROR 1: Unknown color : Menlo-Regular;}"

It seems it is seeing the font type, for some odd reason.


texts: ArcGIS Python API - feature renderer faild in AGOL

texts:
I am trying to save a feature service in AGOL using ArcGIS Python API with symbols. I use the same symbols in a map widget and the features appear correctly but when I save the feature as a feature service in AGOL the symbols are not stored.
I use the following code:
#setting layer symbol properties
layer_options = {
    "type": "FeatureLayer",
    "renderer":"ClassedColorRenderer",
    "title":"{0} Recreation Distribution".format(choice),
    "outline": {
        "color": "White",
        "width": 0.3,
        "type": "esriSLS",
        "style": "esriSLSNull"
        },
    "field_name":"PUD_YR_AVG",
    "opacity":0.6}

#upload shapefile to AGOL
sf_upload = gis.content.add({}, zip_file)
print(sf_upload.id)
published_service = sf_upload.publish()
published_service.update(layer_options)


texts: Error Opening WMS Maplayer: URLError: <urlopen error [Errno 10060]

texts:
I receive the error:

URLError: <urlopen error [Errno 10060] A connection attempt failed because the connected party did not properly respond after a period of time, or established connection failed because connected host has failed to respond>

when attempting to open a webpage containing a WMS map. This is part of a mapping system we run on client machines. I have the system running on multiple clients and have not encountered this error before.

It errors at:

mapServerUrl = settings.MAPSERVERABSOLUTEURL+"?map="+map_file wms = WebMapService(mapServerUrl, version='1.1.1')

It fails a urllib do_open() command, which is not in my code. I tested the map server url manually and received a non-response.

It seems to time-out loading the page. I posted this to the general stack overflow forum, but thought it might be an issue with Mapserver.

The solution from Microsoft for this error is: https://support.microsoft.com/en-us/kb/811087 . This does not work because to the best of my knowledge they are not operating a proxy server. I also tried turning off firewalls, etc because this seemed to work for others having similar issues. Any suggestions for trouble shooting would be helpful.

It has been suggested to increase the timeout, but the server already fails before the default 60 seconds (after about 10 seconds) and I have not experienced this issue on other installations.


texts: Using Python to create new point location based on various azimuth and distances

texts:
I have a feature class with thousands of points. In the same table with these point locations, I have two additional fields. One for azimuth, and one for a length (or distance) to move along that azimuth in order to create a new point location. I am looking for guidance on a python script that would accomplish this task?

table example
x             y          length_fornew_point       azimuth
400460.99 135836.76       5150.04                    55.41
.... ....

Ultimately the script would run through my table, look at the existing coordinates and create a new end point at corresponding specified distance/azimuth in the table.

These two coordinates would then be used to create a new two point polyline.


texts: I have a class that inherits from BaseDynamicLayer and I override the DrawDynamicLayer() method. Is there something special that needs to be done for a successful inDynamicDisplay.DrawPolygon(poly) call? The method crashes every time I call it with a "Exception from HRESULT: 0x80040218." I can use the same point collection to create a Polyline and call DrawPolyline(polyLine) without a problem. I can also successfully use all the other draw commands, like DrawMarker() & DrawRectangle().


texts:
How do I retain symbology that is in QGIS and put it into Mapinfo? I've used the OS Translator II app (https://www.lutraconsulting.co.uk/products/ostranslator-ii/) to apply the styles to the OS dataset in QGIS.

I need to open this dataset in Mapinfo v16 with symbology.

In QGIS i've saved the layer as both .shp and .tab files and neither open with symbology in Mapinfo.

Below is the symbology I want to retain


texts: Extract lines of rice terraces from satellite image

texts:
I have a Pleiades image and multiple rice field boundaries in a shape file.
What I want to do is for every feature in the shape file, clip the Pleiades image and detect the lines of the field terraces. I have already managed to detect the lines using the difference of gaussian blurs, but they look quite messy. Is there a way to do this smarter with better results?
This is the part of my script that detects the lines:
from skimage import filters, exposure

 # I average all bands to a single band
avgs = [ avg_bands(i) for i in clips]

result_rasters = []

# normalize images between 0 and 255
# using scikit-image library for gauss and median filter
for raster in avgs:
    dif_gauss = normalize(filters.difference_of_gaussians(
    raster,
    low_sigma=0,
    high_sigma=3),255)

    med_dg = filters.median(dif_gauss)

    # I look for values smaller than the mean of the image
    med_dg = med_dg < med_dg.mean()

    result_rasters.append(med_dg)


texts: Is it a possibility in OL to add a legend (symbology) directly to (under) particular layer (not only switch on/off the layer)?


texts:
I am trying to compute abs(A-B), where A is a multiple layer raster and B is a single raster.

I have tried the following script into the QGIS Processing Toolbox:

##XXX=name
##Utilities=group
##Raster=multiple raster
##B=raster
##OUT=folder

import glob, qgis
from PyQt4.QtCore import QFileInfo
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer

# Split rasters
layers = Raster.split(';')
output_path = OUT + "/"
suffix = "_suffix.tif"

# Get layer object
layer2 = processing.getObject(B)

for ras in layers:
    # Get layer object
    lyr = processing.getObjectFromUri(ras)
    entries = []

    ras = QgsRasterCalculatorEntry()
    ras.ref = 'lyr@1'
    ras.raster = lyr
    ras.bandNumber = 1
    entries.append(ras)

    raster2=QgsRasterCalculatorEntry()
    raster2.ref='B@1'
    raster2.raster=B
    raster2.bandNumber=1
    entries.append(raster2)

    calc = QgsRasterCalculator( 'abs(lyr@1-B@1)', output_path + lyr.name() + suffix, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
    calc.processCalculation()

for results in glob.glob(output_path + "*" + suffix):
    fileInfo = QFileInfo(results)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

But it doesn't work. When I run the process I get the following message:

unicode cannot be converted to qgis._core.QgsRasterLayer in this context


texts:
I'm trying to plot a bunch of points on a map of Britain and Ireland, and I've been using GeoPandas for that.

I've basically just been following the tutorials to do that, and slicing a world map with a Mercator projection such that it only shows Great Britain and Ireland. But somehow the aspect ratio (projection) ends up distorted.

Here's what my squashed map looks like now:

Here's what I've tried:

Setting the figsize attribute to something with a larger height
Assigning a different CRS
Using a different plotting engine altogether (I can't seem to get Altair to work)

My plot as it stands, with the geopandas code used to generate it, is here on GitHub.

What am I doing wrong, and how do I unsquash this map?