texts: Parallel feature request connection limit

texts:
I have a large pool of feature requests that need to be performed on a PostgreSQL layer. Please assume that these requests cannot be joined into larger requests. If I spawn several threads to perform these requests, there is a limit of 4 database connections to fulfill these requests. All other threads are blocked waiting for one of these 4 connections to open.

Is there a way to increase this limit? Is this limit imposed by QGIS or by the PostgreSQL server?


texts: Having problems installing pyqt on Mac Os sierra

texts:
I am trying to write my first plugin following the tutorial at http://www.qgistutorials.com/en/docs/building_a_python_plugin.html.

Woking on a Mac running sierra.

The problem is that home brew will install only pyqt5 which, so far as I can tell, only works with python 3. At any rate it insists on installing python 3 and then dies trying to install sip but there are no obvious errors in the install log.

I then tried installing pyqt from https://doc.bccnsoft.com/docs/PyQt4/installation.html but again I failed to get sip installed.  The page admits that "installation for mac os can be complex".  I managed to build without problems but running make install tried to install sip in the same directory as python and mac os gave the error "Operation not permitted".  Same thing happens if you try and add any file to that directory.  Yes, I was using sudo ;)  It isn't a straight forward permission problem.

Update:  The fundamental problem is that the python2.7 configure magic does not know about Mac OS SIP (System Integrity Protection) which prevents files in certain directories from being modified even by root.

I have spent a morning hacking the make files produced by various python configure scripts and managed to get sip (not to be confused with SIP, above) installed and to get the pyQT4 Core to find the include files but when I try to make pyQT I get the C compiler complaining about unknown type name which indicates that something is still wrong with the headers.

At this point I give up.  I suspect that the python 2.7 development environment on Mac OS after El Capitan is seriously broken.

I'd be delighted if someone proves me wrong but I suspect the chances of anyone claiming the bonus are near zero.


texts: Drag multiple vertices of different line features in QGIS

texts:
Is it possible to drag multiple vertices of some parallel lines in QGIS at the same time?

I have a line-layer with parallel lines. In some regions I have to move the vertices of all parallel lines so that they still keep their distance to each other.

With the node tool I can only select the vertice of one of my lines I want to drag but not of multiple lines at the same time:

In Multiline-Layers it is possible to drag all lines at the same time but that's not an option because I have to keep all attributes of my single lines.

Are there any other options to drag my lines all at the same time?

Edit1: Just found the CadNodeTool (https://github.com/wonder-sk/CadNodeTool ). With this plugin it is at least already possible to select vertices of multiple lines at the same time. Perhaps this plugin can be extended to move vertices of parallel lines.


texts: After doing TOA correction to a lansat 8 image band using raster calculator, the reflectance is slightly out of range. Why is that so?

texts: the reflectance range was like for example 0.0023 to -1.02356 for some band and -0.0156 to 0.009. Why is this reflectance range slightly out of range? Isn't it should be in exactly in between 0 to 1. If so what can we do to correct this?


texts:
I've been using QGIS 2.14 pretty successfully for about two years. I do also have Arc but tend to use Q more primarily because it did run considerably faster. I then installed QGIS  3 as I noticed it had some really useful functions such as the 'multi edit tool' for the attribute table. However the difference in speed between the two versions is drastic, particularly when working with 'layer properties' or the print composer. More often than not when I open the properties menu my computer freezes and I can go and make a coffee before it's even opened and that's before I've tried to make any changes. Even moving the cursor around at that stage can lead to a 'not responding message'.

As I've not seen any other comments to this effect I can only assume I'm missing something. Perhaps a setting or something. As far as I'm aware I'm not using the program any differently than I was the previous version, but it's currently making a program that transformed the way I work into something that is almost unusable.


texts: Calculated sum of visible areas in print composer

texts: I want to view the sum of area of each habitat displayed in the layout. Only the displayed habitats in the viewport. Which expression do I have to use in the textbox? Now I have calculated them manual. Is it possible to automate this?


texts: How to update a Geoprocessing Service when its Python script changes?

texts:
I want to use a Python script to perform some analysis, and publish this via an ArcGIS Server Geoprocessing Service.

A quick tour of publishing a geoprocessing service says to run the Python script in ArcMap to get a Geoprocessing Result - it's actually this result which you publish as a GP service.

Does this mean that every time you make any changes to the Python script, you need to:

re-run the Python script (in ArcMap) in order to obtain a new Result
publish the Result to ArcGIS Server (over-writing the last GP service)?

I'm hoping there's a faster way to update the script as this will very quickly become tedious during testing...


texts: How to link ArcGIS Online content to my organization's own ArcGIS Online Home Page?

texts:
My organization has an ArcGIS Online account.  We have jumped through all the hoops to set up our Home Page.  However, we now notice that if someone just logs into ArcGIS Online (the main site, not our organization's), while they can search our content, there is no way to automatically link to our ArcGIS Online Home Page.  Any ideas?  We want to direct traffic to our page.

If it is the only option: I have created a hyperlinked thumbnail that leads to our site.  How can I prioritize it so that it always pops to the top of the list when its Tags are searched for?


texts:
I'm trying to make a search engine in OpenLayers 3 by GeoServer features and WFS.

To do this should I use a CQL filter?

This is my connection code if GeoServer mediates OL3 would you be so kind as to show me an example of CQL from OpenLayers.

     vectorWFS = new ol.source.Vector({
            loader: function(extent,resolution,projection){
                $.ajax('http://localhost:8080/geoserver/wfs',{
                    type: 'GET',
                    data:{
                        service:'WFS',
                        version:'1.1.0',
                        request:'GetFeature',
                        typename:'Botanico:Urbana',
                        srsname:'EPSG:3857',
                        outputFormat:'text/javascript'
                    },
                    dataType:'jsonp',
                    jsonpCallback:'callback:loadFeatures',
                    jsonp:'format_options'
                });
            }

        });

        window.loadFeatures=function(response){
            GeoJSON= new ol.format.GeoJSON();
            vectorWFS.addFeatures(GeoJSON.readFeatures(response));
        };

        layerVector = new ol.layer.Vector({
            source: vectorWFS,
            style: new ol.style.Style({
                stroke: new ol.style.Stroke({
                color: 'rgba(0,0,255,1.0)',
                width: 2
                })
            })
            });


texts:
I have the number of travels from one city to nearby cities.
I already created centroids for all of these cities.

And I want to create a map connecting this cities (by lines) and to give those line weight according to number of travels.

How do I create this lines? How do I make the connection between different cities - Table of attributtes?

In summary, I want to end up with something like this map:


texts: Changing object/point to custom symbol based on specific name in MapInfo using MapBasic?

texts:
I am new to programming and Basic.  I only started trying to use MapBasic and have reviewed a few sites and their reference guide to get an idea of what I can do and would like to do.

For this project I am trying to eliminate a monotonous task that is performed each time a user creates a "custom aerial".  Basically we have a table with a range of 15 to 45 (sometimes more and sometimes less) records that are name specific.  For example, let's say they are Fast Food Users (Taco Bell, McDonald's, KFC) Each time we create an aerial we change each record to their respective custom symbol (which I have created) and it uses up a lot of time.  So what I am trying to figure out is how I can write some code that will automate this task by taking my table "retail.tab" and call on a column ("CustImgName") that contains the custom symbol name (TacoBell.bmp, KFC.bmp, etc) and automatically change it to the custom symbol.

In addition, I would also like to find a way to report back if it matches them all or not so that the user can QC their work before proceeding.

I have tried the Hard Code Thematics tool that is out there, but due to the massive amount of records we use in our database it causes performance issues with MapInfo.


texts: Using field calculator to populate one column using other data in same table in ArcMap

texts:
I'm not the best with using the field calculator of ArcMap and I'm unsure of what's wrong with my VBScript

Basically the new column needs to contain data pulled from the same row in the following format
"Name": Streetname-ObjectID

In case it's not clear, the result should be prefaced by the word "Name" and then be populated with the name from the street column, a hyphen, then the corresponding ObjectID.

This is what I've been trying to run more or less, but the VB Script parser won't take it:

'Name': [name] & '_' & [OBJECTID]


texts: Python Error using MMQGIS, global name 'math' not defined

texts:
I am using QGIS 2.16 and trying to compute distances between two points using the MMQGIS plug-in (Hub -distance). This worked like a charm, but recently I have been getting a Python Error whenever I change the measurement unit from "layer unit" to any of the other. The error message is:

"   Traceback (most recent call last):
              File "/Users/henry/.qgis2/python/plugins/mmqgis/mmqgis_dialogs.py",   line 1727, in run
                nameattributename, units, addlines, savename, evenly_distributed, 1)
              File "/Users/henry/.qgis2/python/plugins/mmqgis/mmqgis_library.py", line 4169, in mmqgis_hub_distance
                hubdist = mmqgis_distance(source[0], hub[0])
              File "/Users/henry/.qgis2/python/plugins/mmqgis/mmqgis_library.py", line 544, in mmqgis_distance
                start_lon = start.x() * math.pi / 180
            NameError: global name 'math' is not defined"

Does anybody know how to resolve this issue?


texts: How to Bulk Import to R-tree using a Generator?

texts:
I've been trying to improve one of a couple of scripts I'm writing to be comparable to the other script in performance speed which uses a different approach to tackle the same FGDB problem.  I'm simply trying to find buildings within perennial waters.

It's been brought to my attention even though I'm trying to do an R-Tree via a bounding box spatial analysis, that I may doing it wrong by not doing a bulk import to R-Tree via a generator.

Some concerns I have is it seems I need to manually find each feature's minimum bounding rectangle and also utilize a generator manually defined (which I've copied from the R-Tree website).

My code below which is currently getting me nowhere:

def getMBR(feature):

    global maxX
    global minX
    global maxY
    global minY

    for currFeature in feature[1]['coordinates']:
        if currFeature:
            for x,y,z in currFeature:
                if (maxX == maxY):

                    maxX = x
                    minX = x

                    maxY = y
                    minY = y

                else:

                    maxX = max(maxX, x)
                    minX = min(minX, x)

                    maxY = max(maxY, y)
                    minY = min(minY, y)
    return

def generator_function(someData):

    for i, obj in enumerate(someData):
        yield(i, (obj.xmin, obj.ymin, obj.xmax, obj.ymax), obj)
    return

def applyRTree(layer1, layer2):
    with fiona.open(sourceDatabaseFile, 'r', layer=layer2) as sourceFile2:
        for water in sourceFile2:
            if (('ZI024_HYP' in water['properties']) and (water['properties']['ZI024_HYP'] == 1)):
                waterFeatures.append((water['id'], water['geometry'], water['properties']))
        sourceFile2.close()

    #from rtree import index
    w = index.Property()
    idx = index.Index(properties=w)
    for i, coords in enumerate(waterFeatures):
        getMBR(coords)
        idx.add(i,(minX, minY, maxX, maxY))
        maxX = maxY = minX = minY = float(0.0)

    b = index.Property()
    idx2 = index.Index(properties=b)
    for i, coords in enumerate(buildingFeatures):
        getMBR(coords)
        idx2.add(i,(minX, minY, maxX, maxY))
        maxX = maxY = minX = minY = float(0.0)

    hits = idx.intersection(idx2, objects=True)

try:

        for bldg in buildingLayerList:
            with fiona.open(sourceDatabaseFile, 'r', layer=bldg) as sourceFile1:
                for building in sourceFile1:
                    buildingFeatures.append((building['id'], building['geometry']))
                sourceFile1.close()

            #Use current list of buildings and compare against all hydrography layers
            for hydro in hydrographyLayerList:
                applyRTree(bldg, hydro)


texts:
I have an interpolated temperature data samples that covers a map evenly. The data is in the form of (lat,lon,temp). I want to render that data on top of google maps. I tried google heat maps and the rendering is inaccurate.

I assume one way to do it is to have a layer stored in geoserver. I don't want to use this option. I would prefer if the data are included with the webpage and rendered at runtime.

What technique should I follow?

It seems that is is possible to render polygons in google maps. You could set the opacity, and color. What I did is interpolating temperature at different degrees, and then creating a large number of polygons that get rendered on top of the map.

The color of each polygon is changed based on the approximated temperature value. I got the following results:

Now I have the following isse: File size is too large, about 1MB is needed to store the code of java script.

How can I send the java script code in compact format to the client?

Also let me know if there is an alternative way to render polygons on top of google maps.


texts: Longitudes are out of range when calculating boundry of raster layer

texts:
I want to determine the lat/long coordinates of my raster layer. I calculated them using the geotransform of the dataset.

When I run the script on my file, it produces longitudes that are out of range:

(0.0, 0.0,  '->', 179.88, 90.12)
(0, 751,    '->', 179.88, -90.11999999999998)
(1500, 751, '->', 539.88, -90.11999999999998)
(1500, 0,   '->', 539.88, 90.12)

The longitude ranges from 180 to 540 (it's off by 360)

Here's the code I used:

import gdal, osr

gribdata = gdal.Open('CMC_glb_ABSV_ISBL_250_latlon.24x.24_2018041800_P000.grib2')
gribdata_gtrn = gribdata.GetGeoTransform()

gribdata_bbox_cells = (
    (0., 0.),
    (0, gribdata.RasterYSize),
    (gribdata.RasterXSize, gribdata.RasterYSize),
    (gribdata.RasterXSize, 0),
  )

geo_pts = []
for x, y in gribdata_bbox_cells:
    x2 = gribdata_gtrn[0] + gribdata_gtrn[1] * x + gribdata_gtrn[2] * y
    y2 = gribdata_gtrn[3] + gribdata_gtrn[4] * x + gribdata_gtrn[5] * y
    print (x, y, '->', x2, y2)

The file I'm using is this one, taken from the CMC global forecast.

The latitudes are correct (90 to -90), but howcome the longitude is 180 to 540?


texts: Calculate the coordinate of a point inside a bounding box of a UTM zone map given the bounds, cartesian distance and ellipsoid mesurements

texts:
I am relatively new to the GIS domain.
I recently purchased a 3D model of a city that is mapped to UTM zone 40N, and was provided the coordinates of the 4 corners of the bounding box in Lat Long. In Unity, I am currently trying to calculate the easting & northing of a point on the map (like a mouse click), and vice versa. I am able to get cartesian distance from the corners in Unity meter units, however the resultant coordinate inaccuracy become too large (~30m). I believe I need to use ellipsoid measurements instead, but am unable to figure out how to do so.
Unfortunately all Unity map services are Web Mercator only. I managed to replicate the bounding box in QGIS with UTM zone 40 projection, but couldn't figure out how to port that over into Unity. It is frustrating to see the coordinates I desire displayed on hover, but I have no idea how to retrieve or recreate it in Unity.


texts: I have a trajectory and a map that I have divided in segments. I want to locate each point (using for its id its time-stamp) to the appropriate segment, using qgis. How can I accomplish this?


texts: How I relate Census Blocks to FIPS Division/Subdivisions

texts: I need to be able to relate a Census Block to it's corresponding FIPS Division and Subdivision and I have not found any data file or reference to help me identify that relationship, can you help?


texts: Why are the formats of these coordinates different?

texts:
I'm really sorry if this is a silly question (and I'm sure it is) but I have tried searching but I am a complete newbie to coordinates and maps and I don't know what I'm looking for so am not finding my answer!
I have a gpx file which has coordinates like as follows:

51 47.322,-2 47.218
51 47.340,-2 47.071

I need to add some new coordinates but I am taking these from google maps,but the coordinates are in a different format.
ie:

51.804979,-2.738793

I know there is a logical reason for this as I've read about different coordinate types, but I don't understand why and what types they are and what I need to do to convert from one to the other?
Thanks for any help!
Bex


texts:
I am trying to label counties using choroplethr.  Some posts on the topic exist but do not detail sufficiently how this is done.

My reproducible example/code follows.  How would I change my code so that it displays labels for my counties? If needed, the label variable I'd like to use is "CountyName". Setting c$show_labels to TRUE doesn't seem to work.

library(RCurl)
library(choroplethr)
library(choroplethrMaps)
HSmap<-getURL("https://raw.githubusercontent.com/bac3917/Cauldron/master/NeedsSurv_MAP.csv")
HSmap<-read.csv(text=HSmap)
library(plyr)
HSmap<-plyr::rename(HSmap,c(fipst="region") )             #rename variables
HSmap<-plyr::rename(HSmap,c(ChallSum_mean="value") )      #rename variables

county_choropleth
function (df, title = "", legend = "", num_colors = 5, zoom = NULL)
{
      c = CountyChoropleth$new(df)
      c$title = title
       c$legend = legend
       c$set_num_colors(num_colors)
       c$set_zoom(zoom)
       c$show_labels = TRUE
       Map1=c$render()
}

Map1<-county_choropleth(HSmap,
              title = "2017 Survey - Average Things per Staff, by County",
              legend="Number of Responses",
              state_zoom = "pennsylvania",
              reference_map = FALSE)

Map1


texts:
I'm trying to run the following code to extract different types of climate data for various points:
# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
arcpy.env.workspace = r"X:\AugustComplex"

inRaster = "BurnSev_IA.img"
outPoint = "FirePts.shp"
field="VALUE"

# Execute RasterToPoint
arcpy.RasterToPoint_conversion(inRaster, outPoint, field)

# Bring in the boundary and buffer it by 500m
boundary = "Boundary.shp"
boundaryBuffer = "NewFire_Boundary_500mBuffer.shp"
distanceField = 500
sideType = "FULL"
endType = "ROUND"
dissolveType = "ALL"
dissolveField = ""
arcpy.Buffer_analysis(boundary, boundaryBuffer, distanceField, sideType, endType, dissolveType, dissolveField)

# Clip points to buffer
in_features = "FirePts.shp"
clip_features = "NewFire_Boundary_500mBuffer.shp"
out_feature_class = "FirePts_clp.shp"
xy_tolerance = ""

# Execute Clip
arcpy.Clip_analysis(in_features, clip_features, out_feature_class, xy_tolerance)

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

inPointFeatures = "FirePts_clp.shp"
# Climate variables
inRasterList2 = [["AET.tif", "AET"],
                 ["Snowpack.tif", "Snowpack"],
                 ["CWD.tif", "CWDeficit"],
                 ["PPT.tif", "PPT"]]

# Execute ExtractValuesToPoints
ExtractMultiValuesToPoints(inPointFeatures, inRasterList2)

When I try to run this, I get:
Runtime error  Traceback (most recent call last):
  File "<string>", line 7, in <module>
  File "e:\arcgis\desktop10.7\arcpy\arcpy\sa\Functions.py",
    line 7407, in ExtractMultiValuesToPoints     bilinear_interpolate_values)
  File "e:\arcgis\desktop10.7\arcpy\arcpy\sa\Utils.py",
    line 53, in swapper     result = wrapper(*args, **kwargs)
  File "e:\arcgis\desktop10.7\arcpy\arcpy\sa\Functions.py",
    line 7402, in Wrapper     bilinear_interpolate_values)
  File "e:\arcgis\desktop10.7\arcpy\arcpy\geoprocessing\_base.py",
    line 510, in <lambda>     return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 010555: Internal error in Sample: %s.
Failed to execute (ExtractMultiValuesToPoints).

It might be worth noting that the block before that uses the same function and works fine. I haven't been able to find much about what this error means.
Any suggestions for what to try?


texts: Creating new highway mapping data? Help!

texts:
I'm hoping somebody may be able to help or at least point me in the right direction, with a GIS/Mapping query that is a little unusual.

I have been tasked at work with looking at how we create a Google Maps style system for our internal roadway system - this roadway covers around 45 miles and the level of detail exceeds that supplied by any available mapping data (OS, Google Maps, OSM) source.  I've looked at most the available data and its sources and its just not an option - therefore we need to have maps drawn/created to our specific requirements, which may seem OTT, but its the only way. Also this is not detail that can be done as an overlay on existing map data, the 'bosses' have already ruled that option out.

In the main that is easily done in a variety of software.  However, the company want it available, link to GPS data for junctions, signposts etc, and available in an online Google Maps style interface, which means that once created, the map needs to be rendered into tiles. My first thoughts was a combination of CloudMade and the Leafletjs javascript library, but the former does not appear to allow you to import your own mapping data.

I've searched extensively on the net and can find numerous options for overlaying additional information on existing data, but nothing on how maps can be drawn and rendered, then tiled and displayed.  Leafletjs solves the latter issue re-displaying, but can anyone point me to software that can create maps, link specific points to GPS co-ordinates, and export the whole thing through tiles etc..?

Are maps drawn for this purpose, or do you start with a GPS 'track' of the roadway and covert it/update/render it to a map?

I'd welcome any thoughts/suggestions as to how the project could proceed as its not my field of expertise. Help!


texts:
for personal usage, I downloaded an openstreetmap map (as a *.osm) and try to modify it using  the Josm editor.
I want to trace a path in the card I walked and want to remove unnecessary things from this map. Other paths are highlighted, because they are access-restricted. This is saved in the preset Highways/Streets/Road Restrictions.

 This is how Josm displays what I mean.

How can I delete this preset or hide the highlighting? I don't want to delete the whole path, though.

For better illustration, path I mean goes from the upper-left corner to the right and is highlighted by the application in blue.


texts: Why do my XY cordinates not make sense / show up in ArcGIS 9.3?

texts:
In the Add XY Data menu I selected to add a table containing UTME and UTMN coordinates. I selected a WGS84 projection b/c that's what the GPS device that meade the readings was set to use. However, my raster map (and other features) in ArcGIS use a NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet projection. As a result (I assume), when I add my XY data layer, it doesn't show up anywhere on the existing map.

My first attempt to fix this was to re-project the new XY points:

I tried the Project tool, and chose to project from WGS84_UTM_Zone_17S to NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet using NAD_1983_To_WGS_1984_5 as my geographic transformation.
The values of my X and Y columns did not change, and the points still don't seem to exist on my map.

My second attempt to fix this was to convert my new XY points from meters into feet (US) [since the other layers are in a 'feet' pojection.

I tried converting this reprojected layer from meters to feet (Using the "Calculate Geometry" option in the attribute table).
It converted my UTME to a sensible number (from 689583.7 to 2744074), but the UTMN conversion makes no sense (from 4008525 to -52330979).

What am I doing wrong and/or what can I do to make this work?

Note:nearby points already on my existing map have coordinates {2032037.05127, 892767.424497}

Below are the properties of (first) my added XY points, (second) my underlying raster DEM map, and (third) the foot-converted reprojected XY data.

Data Type: XY Event Source
Location: C:\Users\paynec\Desktop
Table: H23 Geocoordinates2.csv
X Field: UTME
Y Field: UTMN
Has Object-ID Field: No

Projected Coordinate System:    WGS_1984_UTM_Zone_17S
Projection:          Transverse_Mercator
False_Easting:       500000.00000000
False_Northing:      10000000.00000000
Central_Meridian:    -81.00000000
Scale_Factor:        0.99960000
Latitude_Of_Origin:  0.00000000
Linear Unit:         Meter

Geographic Coordinate System:   GCS_WGS_1984
Datum:          D_WGS_1984
Prime Meridian: Greenwich
Angular Unit:   Degree

###############################################

Spatial_Reference:    NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
Linear_Unit:          Foot_US (0.304801)
Angular_Unit:         Degree (0.017453292519943295)
False_Easting:        2000000.002616666
False_Northing:       0
Central_Meridian:     -79
Standard_Parallel_1:  34.33333333333334
Standard_Parallel_2:  36.16666666666666
Latitude_Of_Origin:   33.75
Datum:                D_North_American_1983

###############################################

Projected Coordinate System:    NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
Projection: Lambert_Conformal_Conic
False_Easting:  2000000.00261667
False_Northing: 0.00000000
Central_Meridian:   -79.00000000
Standard_Parallel_1:    34.33333333
Standard_Parallel_2:    36.16666667
Latitude_Of_Origin: 33.75000000
Linear Unit:    Foot_US

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: Merge overlapping/neighboring polygons - GeoPandas

texts:
I am reading a shapefile and I have extracted a specific row from it that is a Multipolygon. I can extract all the polygons in it with .explode() which works great and I get 130+ polygons.

How can I merge the ones that are next to each other?

I know the 130+ resulting polygons represent two islands, and I want to get each as a polygon. I have looked up how to do it but haven't been able to find anything and this is my first project ever using geopandas and shapely.

gdf=gpd.read_file("comunas.shp").to_crs({'init': 'epsg:4326'})
Rapa_Nui=gdf.loc[gdf["Comuna"]=="Isla de Pascua"]
all_Rapa_Nui=Rapa_Nui.explode()


texts: Converting Full Motion Video for ArcGIS Desktop

texts:
I am experimenting with Full Motion Video in ArcGIS Desktop.

I used tutorial data to learn more about Full Motion Video, and then acquired UAS video from my employer.  These videos did not work, and after investigation I found that this was due to there being 12 audio streams in the videos.  After getting rid of the audio using a video conversion application to get rid of the audio, I could sort of use it.  In ArcMap I could see the video with the FMV plugin, but then realized there was no metadata in the converted video.

I would like to know how I could convert video using COTS software, while keeping the metadata, but losing the audio, and how I could extract the metadata from the original video to save it as csv, and then recreate the video using multiplexer to combine the video and csv.

I am using ArcGIS Desktop, FMV version 1.3.2


texts: Choose Layer by name and calculate all polygons of this layer by adding a new column "Area"

texts:
I'm new to QGIS AND Python.

This is how I'm trying to calculate all areas of the layer named polygons:

def add_area():

layers = QgsMapLayerRegistry.instance().mapLayers().values()
foundlayer = False
for layer in layers:
    if layer.name() == 'polygons':
        inputlayer = QgsMapLayerRegistry.instance().mapLayersByName("polygons")[0]
        print ('polygons-layer found')
        iface.setActiveLayer(inputlayer)
        foundlayer = True
if not foundlayer:
    print ("polygons does not exist")
    #maybe not the most elegant way, but it works till here as it should :-)
area = 0
#this creates me a new column named "Area"
provider = inputlayer.dataProvider()
provider.addAttributes([QgsField("Area",QVariant.Double)])
#and now I'm iterating through all features to calculate each polygons area...
features = inputlayer.getFeatures()
f = features.next()
for f in features:
    attribute = (f.attributes())
    print(attribute)
    calculator = QgsDistanceArea()
    area = calculator.measurePolygon(f.geometry(). asPolygon()[0])
    calculator.convertMeasurement(area, QGis.Degrees, QGis.Meters, True)
inputlayer.updateFields()

Even it creates the new column in the attribute table, I'm getting this Error:

Layer polygons: SQLite error: duplicate column name: Area SQL: ALTER TABLE "polygons" ADD COLUMN "Area"

tried to google it but couldn't find a answer in any forum, that I understand and that help's me solving it.

Do I get the message, because the column "Area" is already created? (that is possible, because I naturally tried my code several times! But I got this message already the first time...!?)


texts:
I'm trying to add multiple images to OpenLayers. I have the following (some code is excluded):

var images = [image1,image2,image3] // these are url images on my server
var image = new OpenLayers.Map($('#image'),
               {
                   units: 'pixels'
               });
// add a bunch of controls including layer switcher.
var vectorLayer = new OpenLayers.Layer.Vector("Vector Layer");

for(var i =0; i < images.length; I++)
{
  var info = getInfo(images[i]);
  var viewPort = image.getViewport();
  var viewHeight = jQuery(viewPort).height();
  var viewWidth = jQuery(viewPort).width();
  var graphic = new OpenLayers.Layer.Image(
                        info.name,
                        images[i],
                        new OpenLayers.Bounds(0, 0, info.width, info.height),
                        new OpenLayers.Size(vw * 0.9, vh * 0.9)
                    );
  image.addLayers([graphic,vectorLayer]);
  image.zoomToMaxExtent();
}

This is a fairly simple example of what I am doing and I can see the images are getting loaded. All three images are different sizes (small, medium, large).

Once everything is loaded and I can see the last image loaded perfectly:

I switch to the middle image the scale and positioning are all wrong.
I switch back to the last image and the scale is right.
I zoom into the image (last one)
I switch back to the middle image and I can't see anything any more.

I'm not too sure if what I want to do is even possible, but here goes. I want to load all images but have them all independent. So when I zoom into one image the state isn't the same on another.

Please can someone offer advice?


texts: How to pre-configure GeoServer in docker image?

texts:
I am using GeoServer 2.16.0 docker image build by myself based on kartoza's dockerfile.
The only pre-configuration is done via environment files, e.g. geoserver.env:
ENABLE_JSONP=true
FOOTPRINTS_DATA_DIR=/opt/footprints_dir
GEOSERVER_ADMIN_PASSWORD=myawesomegeoserver
GEOSERVER_DATA_DIR=/opt/geoserver/data_dir
GEOSERVER_NODE_OPTS="id:$hostname - ($host_ip);background:red;color:white"
GEOWEBCACHE_CACHE_DIR=/opt/geoserver/data_dir/gwc
#JAVA_OPTS="-Duser.country=US -Duser.language=en"
MAX_FILTER_RULES=20
OPTIMIZE_LINE_WIDTH=false

I would like to configure more, e.g. workspaces, stores, but cannot find any according documentation.
Is there any better way than the following?

Configuring a running container having the DATA_DIR mounted
Create an own image based in kartoza's and add the DATA_DIR content?

Solution (16.02.2021)
I added the following two lines to my dockerfile:
COPY ./data_dir_init /tmp/data_dir_init
COPY ./init-workspace.sh /scripts

The init-workspace.sh contains the following:
SETUP_LOCKFILE="${GEOSERVER_DATA_DIR}/.init_data_dir.lock"
if [ -f "${SETUP_LOCKFILE}" ]; then
    echo -e "\n\n Skip GeoServer init - already done\n\n"
    exit 0
fi

if [ "${DEBUG}" ]; then
    set -e
    set -x
fi

echo -e "\n\nStarting init of GeoServer data dir\n\n"

cp --verbose --recursive /tmp/data_dir_init/{styles,workspaces} $GEOSERVER_DATA_DIR
cp --verbose /tmp/data_dir_init/{csw,global,logging,wms,wps}.xml $GEOSERVER_DATA_DIR
mkdir--verbose --parents $GEOSERVER_DATA_DIR/data

if [ ! ${DEBUG} ]; then
    rm --recursive --verbose /tmp/data_dir_init
fi

# Put lock file to make sure data directory is not reinitialized on restart
touch ${SETUP_LOCKFILE}

The folder data_dir_init contains all content from the first session, when I startet an empty geoserver instance and created the base set-up I needed.


texts:
I've got a lot of points with Z-values, where I would like to construct their extent using IDW (or similar). I would like to see their extent inside my buffer, is this possible or is the tool not designed to move beyond/behind the last/first point? See picture for explanation, I can use the buffer to cut the one in shown when running the algorithm, but can't get the extent quite right.

The points are spaced approx 10 m from each other, the buffer is 50 meters from centerline.

I've got access to 3d analyst, spatial analyst and arcgis advanced desktop.


texts: Changing legend label when cycling through Atlas created with QGIS composer

texts:
I have project which gets updated regularly and can have multiple areas which need to be mapped individually to go into separate reports. I generate an atlas in the composer and I have created a rule-based style using the following:

attribute(  $currentfeature , 'id') =  attribute(  @atlas_feature , 'id')

This allow the user to go through each property, and if they are close together, then you only see the features for that specific area, not adjacent ones that might be in the map view. In the map composer, I can get the name of the area I am mapping and use in the title of the map using the string below (which changes each time you go to the next area within the atlas.

[% @atlas_pagename %] Rest of Title

One of the key issues i'm facing now is showing in the legend the areal measurement of what I am mapping. The following code; if I use a categorized style; will give me that label I am looking for (Property Boundary 100ha), but as it is not rule based, the surrounding features are not filtered out based on the atlas.

  concat('Property Boundary ',round(area(transform($geometry,'EPSG:4283','EPSG:3308') )/10000), 'ha')

I tried using a rule based system and then applied a categorised style, but the label came out as the following in the label.

concat('Property Boundary ',round(area(transform($geometry,'EPSG:4283','EPSG:3308') )/10000), 'ha') = 'Property Boundary = 100ha'

Is there a way I can filter polygons based upon the atlas and dynamically change the label/legend value to reflect the area using the tools within QGIS (without running Python scripts).


texts: Convert GeoJSON labels MBR to MapInfo text object

texts:
I'm using Dymo (https://github.com/migurski/Dymo/) to generate automatic point labels placement. The output of Dymo is GeoJSON containing MBRs of the positioned labels with the names of the labels in the attributes table. I would like to convert the GeoJSON table to a MapInfo text object table to finish the placement manually in MapInfo.

I can't find any info about that question.

Edit: Sample of Dymo output in JSON:

{
  "type": "FeatureCollection",
  "features": [
{
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -74.033500000000004,
          45.677299999999995
        ],
        [
          -74.028294368743914,
          45.677299999999995
        ],
        [
          -74.028294368743914,
          45.677506895418496
        ],
        [
          -74.033500000000004,
          45.677506895418496
        ],
        [
          -74.033500000000004,
          45.677299999999995
        ]
      ]
    ]
  },
  "type": "Feature",
  "properties": {
    "preferred_placement": "",
    "overlaps": 0,
    "font_size": 6,
    "name": "Label Name",
    "font_file": "Fonts/arialbd.ttf",
    "point_size": 0,
    "justified": "left"
  }
},
{
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -74.030500000000004,
          45.768493442022979
        ],
        [
          -74.025878017425526,
          45.768493442022979
        ],
        [
          -74.025878017425526,
          45.768700000000003
        ],
        [
          -74.030500000000004,
          45.768700000000003
        ],
        [
          -74.030500000000004,
          45.768493442022979
        ]
      ]
    ]
  },


texts: Varying page size of maps within an Atlas in QGIS?

texts:
I would like to create an atlas of a series of maps covering a number of woodlands and to have them all at the same map scale.
This would work best visually if some of the maps were in landscape, and some were in portrait orientation.
Is it possible to do this, or do all maps within an atlas have to have the same paper size and orientation?


texts: Change contour line stroke weight based on elevation?

texts: I have shapefiles of contour lines that I have imported in Qgis. At the moment the contour lines are all the same stroke weight. I'm a beginner to this program, but wanted to know if its possible to change the stroke weights of the contour lines based on the elevation?


texts: r.walk and r.external grass function problems in python console

texts:
I would like to run the r.walk grass function with these codes:

import processing
from qgis.core import QgsRasterLayer
from PyQt4.QtCore import QFileInfo
#processing.alghelp("grass:r.walk")
path= 'd:/Works/_EvacuationSimulation/_raster/walk/'
filename1 = 'elevation.tif'
filename2 = 'friction.tif'
filename3 = 'destination.shp'
filename4 = 'walked.tif'

elevation = QgsRasterLayer(path + filename1, 'elevation')
friction = QgsRasterLayer(path + filename2, 'friction')
destination = QgsVectorLayer(path + filename3, 'destination',"ogr")

extent = elevation.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

min_area = 0.000100
#pixelsize = 0.000000
pixelsize = 0.100000

#processing.runalg('grass:r.walk', elevation, friction, destination, "0", "100", "", "0.72,6.0,1.9998,-1.9998", "1.0", "-0.2125",True,True, True, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), pixelsize, "-1", min_area, path + filename4)
processing.runalg('grass:r.walk', elevation, friction, destination, "0", "100", "", "0.72,6.0,1.9998,-1.9998", "1.0", "-0.2125",True,True, False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), pixelsize, "-1", min_area, path + filename4)

After an execution without getting a result raster, the console log shows problematic rows:

    C:\PROGRA~1\QGIS2~1.14\bin>r.external -r input="d:\Works\_EvacuationSimulation\_raster\walk\unit-elevation.tif" band=1 output=tmp1482592256518 --overwrite -o

    ERROR: r.external: Sorry, <r> is not a valid flag

    C:\PROGRA~1\QGIS2~1.14\bin>r.external -r input="d:\Works\_EvacuationSimulation\_raster\walk\friction-cost.tif" band=1 output=tmp1482592256519 --overwrite -o

    ERROR: r.external: Sorry, <r> is not a valid flag

    C:\PROGRA~1\QGIS2~1.14\bin>r.walk  elevation="tmp1482592256518" friction="tmp1482592256519" start_points="tmp14825922565110" max_cost="0" percent_memory="100" walk_coeff="0.72,6.0,1.9998,-1.9998" lambda="1.0" slope_factor="-0.2125" -k -n output=outpute943ea890dd94b879185c7266a98f4c3 --overwrite
    GRASS_INFO_ERROR(4668,5): Raster map <tmp1482592256518> not found

    GRASS_INFO_ERROR(5548,1): Raster map <outpute943ea890dd94b879185c7266a98f4c3> not found

Probably the problem comes from the built-in r.external function of the r.walk script, which uses a no longer exist '-r' flag.

Can you help to solve this?


texts: Raster Time Series, express years in time series

texts:
What's the best option to express years in Raster Time Series object in R.
I define years as dates:

YEARSSEQ <-seq(as.Date("2008-01-01"), as.Date("2009-01-01"), "year")

and then I pass them in rts constructor function

myrts<-rts(myrasterstack,YEARSSEQ)

Passing  as

YEARSSEQ = c("2008", "2009")

returns the error

  time must be one of time-based classes including: 'POSIXct', 'POSIXt', 'Date', 'timeDate', 'yearmon','yearqtr'

Is there any better option?


texts: Suggestions on how I can depict numerous cell towers with different ranges?

texts:
I am trying to make buffers dependent on cell tower ranges. Does anyone have any suggestion on how to insert distance buffers of various that are dependent on attribute values in one layer?

I am using QGIS 2.18 on a MacBook

Am I understanding correctly: "various distance buffer" can make different buffer sizes dependant on feature values?

This is what i did:

  selected tower layer > vector > geoprocessing tools > variable
  distance buffer > input layer: towers > distance field:
  3499_tower_range_mi > segments: 1609.344 because value in meters is
  about 1 mile > Run

That did not work, I got a message saying "float() argument must be a string or a number"

I am not sure what that means. My attribute column for the tower's ranges is "double | Real"

How can I change this to be an integer value?

Or does anyone have any suggestions on how I can depict numerous cell towers with different ranges?


texts: Where to store line events associated with Pipeline in UPDM(2019)?

texts:
I am migrating pipeline data to UPDM(2019). The current pipeline database has geometric network that needs to be migrated to Utility Network in UPDM. I have mapped the relevant feature classes i.e pipeline feature class containing the pipeline' geometries and attributes is mapped to PipelineLine feature class in UPDM. But there are other tables in my database that have the linear line events on those pipelines.
As per the document "it is possible to replace the Centerline featureclass with the PipelineLine featureclass", so my question is if the PipelineLine feature class has the actual geometry of the pipelines where would I store line events in UPDM?


texts: Derive elevation from rgb raster using sampled points

texts:
I have an RGB elevation raster from a seabed. This raster is interpolated from points which I also have access to. I don't want to download all points and do an interpolation myself. Now the problem is the interpolated raster is only containing RGB values and no elevation. I want to add those values in a new band using (hopefully) only a small portion of the points. In other words; can i derive the elevation from the rgb values using some samples easily? Preferebly in qgis and maybe some python.
Why? So i can change the colours to my wishes and maybe visualize it in 3D


texts: How to activate snapping for layer you are drawing polygons in?

texts: I want to draw multiple polygons in a layer that have boundaries without a gap, so I activated snapping for that layer. But it doesn't work. For other layers it works fine.


texts: Repeating steps based on different field using ModelBuilder?

texts:
I have built a model to interpolate air pollution from 40 recording stations using IDW.

The results are then exported to an average for census divisions.

I need to repeat this process hundreds of times for every month for 15 years for five different types of air pollution.

I have all of the data by year as a new field in the recording stations' file.

Table:

Recording_Station 03jan 03feb 03mar 03apr 03may 03jun 03jul ....
1                  26.5 27.5  22.2  32.3  34.5  45.4  55.4
2                  19.1 12.2  45.3  33.3  36.6  55.4  65.5
3                  11.1 15.6  44.4  34.4  66.6  55.4  44.6

etc.

How can I re-run the model for every field without manually changing them each time?


texts: Merging raster files produces error in QGIS?

texts:
I'm trying to merge (Raster -> Miscellanous -> Merge) four raster files into one. Each file has same CRS and are saved as .tif files. Merge produces error:

"The process crashed some time after starting successfully."

The layers themselves don't overlap.


texts:
I need to make skyline for my DEM with buildings on it. The height of building is height of DEM + height of building.

The problem that tool select me improbable zone:

It's not possible that part of buildings (polygon have fixed height) are visiable and part of not.

ArcGIS 10.4.1
Radius 100 meters

I am using skyline tool.


texts: Calculating intersect of two shapefiles using QGIS

texts:
I have two shapefiles.
1. the clinic buffers 500m area.
2. the postcode area.

Is there a simple way to calculate the intersect two shapefiles ?
I'm using QGIS Dufour. I tried to use fTools to intersect. However, it didn't show me actually how much every piece of the green area. It just combined the fields of two shapefiles.

How do I calculate the intersect of two shapefiles?

I used "Field Calculator" to calculate the intersection with $area. But those numbers are bigger than original area.


texts: Feature-level authorisation in ArcGIS Server REST services

texts:
Is feature-level authorisation even remotely possible when accessing ArcGIS (10) via the REST APIs? From what I understand I can secure individual services but at most this gives me control over who can read / write data across the service, and doesn't cover individual layers or features. Is this the case?

I was surprised that there are lots of REST-based (e.g. JavaScript) examples on editing but nothing obvious that looks at who is allowed to edit data - a primary concern in any real-world application.

If feature-level security is not possible what are my options? So far it seems like I must implement my own editing service (possibly through a Server Object Extension?) which can consult Access Control Lists whenever an edit operation is requested. I imagine I could reuse a lot of what's already in the JS API for editing via the Geometry Service but just point it at my own service and fill in the gaps.


texts: Triggering built-in events programatically OpenLayers

texts:
Is there a way to trigger openlayers event from console. In this case:
 ol.interaction.Select event which listens 'select' (ol.SelectEvent.select). It triggers when feature is selected or deselected by mouse click.


texts: Getting my expression to process the text in field

texts:
I am building an expression, whose function is to return True if a characteristic is selected.
The expression requires two arguments, the Feature and the name of the layer.
If I enter the name of the layer as text it works properly. But if I enter a field containing the name of the layer it does not work, even if the record contains exactly the name.
Here is the code of the expression:
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args='auto', group='Custom')
def isSelected(entidad, capa, feature, parent):
    proy= QgsProject.instance()
    capap= proy.mapLayersByName(capa)[0]
    les= list(capap.getSelectedFeatures())
    boolean=False
    for f in les:
        if f.id() == entidad.id():
            boolean= True
            break
        else:
            continue
    return boolean

Note:

The equivalent function, exists (is_selected) but does not work properly in my QGIS, version 3.4
The code works properly on the console


texts: Rasterize a shapefile with Geopandas or fiona - python

texts:
I need to rasterize a really simple shapefile a bit like this http://tinyurl.com/odfbanu. Which is just a shapefile coutaining counties in the US.I have seen this previous answer : GDAL RasterizeLayer does not Burn All Polygons to Raster? but I was wondering if there is a way to do it using Geopandas or fiona and maybe rasterio for the tiff writing part.

So my goal is to rasterize it and assign a value to every polygon sharing a comon value, LSAD in the exempe.

So I wrote the beginning of the code inspired by shongololo in the thread :
Dissolving polygons based on attributes with Python (shapely, fiona)?.

from geopandas import GeoDataFrame

name_in = 'cb_2013_us_county_20m.shp'

#Open the file with geopandas
counties = GeoDataFrame.from_file(name_in)

#Add a column to the Geodataframe containing the new value
for i in range (len(counties)):
    LSAD = counties.at[i,'LSAD']
    if LSAD == 00 :
        counties['LSAD_NUM'] == 'A'
    elif LSAD == 03 :
        counties['LSAD_NUM'] == 'B'
    elif LSAD == 04 :
        counties['LSAD_NUM'] == 'C'
    elif LSAD == 05 :
        counties['LSAD_NUM'] == 'D'
    elif LSAD == 06 :
        counties['LSAD_NUM'] == 'E'
    elif LSAD == 13 :
        counties['LSAD_NUM'] == 'F'
    elif LSAD == 15 :
        counties['LSAD_NUM'] == 'G'
    elif LSAD == 25 :
        counties['LSAD_NUM'] == 'I'
    else :
        counties['LSAD_NUM'] == 'NA'

Really easy stuff so now I am wondering how can I actually write those shapes to a tiff. I began working with Geopandas as I believed that was the best option but if you have an fiona suggestion I am up for it too.

I found a piece of code from rasterio which seems to be able to take a shapely geometry and burn it into a new raster http://tinyurl.com/op49uek

# I guess that my goal should be to load my list of geometries under geometry to be able to pass it to rasterio later on
geometry = {'type':'Polygon','coordinates':[[(2,2),(2,4.25),(4.25,4.25),(4.25,2),(2,2)]]}

with rasterio.drivers():
    result = rasterize([geometry], out_shape=(rows, cols))
    with rasterio.open(
            "test.tif",
            'w',
            driver='GTiff',
            width=cols,
            height=rows,
            count=1,
            dtype=numpy.uint8,
            nodata=0,
            transform=IDENTITY,
            crs={'init': "EPSG:4326"}) as out:
                 out.write_band(1, result.astype(numpy.uint8))


texts: How to set Merge output as Select Input in ModelBuilder?

texts:
I'm using ModelBuilder. I tried to take the Output from the Merge tool as Input to the Select tool. But it doesn't work.
I'm trying to Merge 90 Shapefile lines into one. Then "Select" some features (with an expression) from the WHOLE one to do more things.


texts: ArcMap tool for deleting polygon with C#

texts: I need to create custom tool for ArcMap, which would let the user delete polygons by mouse clicks. I am pretty new in this thing, so how can I retrieve a selected polygon so I could delete it? Plus do I need to load some sort of function to renew the view?


texts: Labels disappear zooming IN in Qgis 2.2 Composer?

texts:
I have a simple map with a polygon Layer of 2 Polygons. Labeling is turned on, so is the option of "Show all labels for this layer, even if they collide".

In the print composer the labels keep disappearing if i zoom in! Also when I export as image (I tried with .bmp, jpeg and pdf) only one of the two polygon label shows up!?!

What the hell is going on? I get that sometimes labels dissapear when u reduce the scale and the labels run out of space, but by zooming in is new to me!

Is there an option in Qgis or the Qgis composer that I am overseeing?

Please help, this is driving me nuts!

Regards,
Martin


texts: gdal_viewshed's results are suffering some glitches vs visibility analysis from SAGA

texts:
I'm seeing some artifacts when using gdal_viewshed vs the ta_lighting tool from SAGA-GIS. Small differences are noticeable between both visibility results.
As I wasn't sure which one was "wrong" relative to the other, I decided to trace a profile with QGIS on both visibility rasters. The observer position and the DEM used in both computations are exactly the same. The DEM is in a Cartesian space.
The SAGA result seems much more accurate as, in between the two grayed blocks which represent shadows in the SAGA results on the figure hereunder, gdal is telling me that this region, which should clearly be visible, is not. Both runs are what I would call "vanilla computations", i.e. no fancy options other than the mandatory ones were used.

The highest terrain peak is at only ~580 meters from the observer on this figure. And as you don't get fooled by the vertical exaggeration in the graph, the region is relatively flat, which I guess may somehow also be an issue. I also have the feeling that gdal is doing some simplifications somewhere, because it's really much faster than SAGA's computation (less than 5 minutes vs more than an 1 hour on the same machine).
Does anybody know, by any chance, what may cause these differences? Especially what may be wrong in the GDAL based tool which I used this way:
$ gdal_viewshed -f "GTiff" -oz 5 -ox ${X} -oy ${Y} -md 10000 -vv 1 -ov -99 input_dem.tif visibility.tif

(X and Y are both integer values)
I also thought that it's because gdal is using a default coefficient to consider the effect of the curvature and refraction and not SAGA. So I ran the gdal computation again using -cc 0 to discard this; absolutely no differences from the run without using this option. (Just in case, I also run it using -cc 1; this actually induces some differences, but not along the view-sight line depicted on the above graph.)

Source: https://gdal.org/programs/gdal_viewshed.html#cmdoption-gdal_viewshed-cc

SAGA version is 7.1.0
GDAL version is 3.3.0dev from the Docker image osgeo/gdal:ubuntu-small-latest
OS is Ubuntu 18.04.


texts:
I am using R to work with multiband raster files and am running into a problem. I have named the bands, but when I export in netcdf format, the band names are being stripped. Anyone know how to specify band names when exporting netcdf files?

For example, if I have a raster object named "ten", and I type names(ten) I get the vector of band names. If I write ten to disk and then reimport using

newraster = brick("ten.nc"), then the names are stripped away. However, brick() is able to import other .nc files fine with the names attached.

The code I am using is

z = writeRaster(ten,filename = a,overwrite=TRUE,) ### Output to ./outdir folder. This is monthly data


texts: Can GeoServer GetFeatureInfo return class name instead of raster pixel value?

texts: I have a land use / cover raster map published in GeoServer 2.6. When I use the query tool I get the pixel value (eg.: 1). Is there a way for geoserver to return what the value 1 means in my land use raster? For instance, 1:Agriculture; 2:Forest etc...


texts: Convert string column back to geometry column?

texts:
I had a table with a column set to 'point' in doctrine2, but this was changed to 'string' and updated.

I can't change it back to 'point' in doctrine2, I get an error:

[Doctrine\DBAL\Exception\DriverException]
  An exception occurred while executing 'ALTER TABLE opr_records ALTER the_geom TYPE public.geometry':
  SQLSTATE[42804]: Datatype mismatch: 7 ERROR:  column "the_geom" cannot be cast automatically to type geometry
  HINT:  Specify a USING expression to perform the conversion.

[Doctrine\DBAL\Driver\PDOException]
  SQLSTATE[42804]: Datatype mismatch: 7 ERROR:  column "the_geom" cannot be cast automatically to type geometry
  HINT:  Specify a USING expression to perform the conversion.

[PDOException]
  SQLSTATE[42804]: Datatype mismatch: 7 ERROR:  column "the_geom" cannot be cast automatically to type geometry
  HINT:  Specify a USING expression to perform the conversion.

There is already data in the column for example:

id,geom
0,0101000020E610000091AE217D7A881A40551A3D49BD424740
1,0101000020E61000009F6A8C15600F1B407EDF9D1DCA3E4740
2,0101000020E6100000525C9D4D06831A40F459A7CCFF424740

I could connect to the DB and perform the update manually, however I am unsure on the SQL to perform ? Something like:

ALTER TABLE opr_records ALTER COLUMN the_geom TYPE geometry USING ST_geom(the_geom);


texts:
I closed the layers panel and don't know how to open it back up.
Could someone give me a simple explanation as to which tab it is located on so I can open it back up?


texts: Export to CAD from Arcgis while preserving the symbology

texts: Is there a way or tool that you can export layers (feature classes - points, lines and polygons) from Arcgis desktop 10.1 into a DWG/DXF, while still preserving the symbol (basic symbol - colours and basic shapes) from Arcgis?


texts: ArcGIS project WGS_1984_UTM_Zone_32N to GCS_WGS_1984

texts:
I'm new to GIS.

I have a shapefile for which I defined the coordinate system (using the Define Projection tool) to WGS_1984_UTM_Zone_32N and now
I'm trying to reproject the shapefile from WGS_1984_UTM_Zone_32N to GCS_WGS_1984.

I'm using the Project tool but I always get the following error and the projection never happens.

ERROR 999999: Error executing function
invalid extent for output coordinate system.
The operation was attempted on an empty geometry

I run the Repair Geometry tool ,it works fine, but then when I try to reproject to GCS_WGS_1984 I still get the same error as above

I don't understand what I'm doing wrong.

I know similar questions have been asked in the past but I couldn't find one where someone gets an error.


texts: GeoServer/WMS: get original raster image in highest resolution

texts:
I am using the WMS interface of GeoServer to retrieve raster data. The raw raster data is stored in an Oracle database,and based on the WMS BBOX parameter and height/width GeoServer decides which zoom level and resolution/pyramid level will be used - as far as I understand.

For my requirement I need to define a BBOX and as a result I want an image (e.g. GeoTIFF) with the highest possible resolution (like the resolution of the raw images in the database) returned. I don't care about the file size, height and width of the image.

Are there WMS parameters to accomplish that?


texts: GeoServer styles AJAX requests are blocked due to X-FRAME-OPTIONS

texts:
I am running GeoServer behind a Nginx proxy with an SSL cert. Everything is working well (that I can tell so far) with the exception of editing layer styles.

Performing any action that involves AJAX requests on that page results in the requests being blocked due to X-FRAME-OPTIONS.

Refused to display 'https://myserver.com/geoserver/web/wicket/bookmarkable/org.geoserver.wms.web.data.StyleEditPage?10-1.IBehaviorListener.0-styleForm-context-tabs~container-tabs-1-link&name=State+Boundaries&wicket-ajax=true&wicket-ajax-baseurl=wicket%2Fbookmarkable%2Forg.geoserver.wms.web.data.StyleEditPage%3F10%26amp%3Bname%3DState%2BBoundaries'
in a frame because it set 'X-Frame-Options' to 'DENY'.

I have disabled all X-FRAME-OPTIONS that I could find, and from what I was reading, GeoServer should default to SAMEORIGIN. Any attempts I made at overriding the policy in /webapps/geoserver/WEB-INF/web.xml prevented GeoServer from even starting up.

If I add my X-FRAME-OPTIONS declaration in Nginx back in, it also gets denies, saying there were multiple options set, an it fell back to DENY.

I'm quite certain this has to do with the Nginx proxy (since GeoServer runs just fine without it):

location / {
    proxy_pass http://127.0.0.1:8080/geoserver/;
    proxy_pass_header Set-Cookie;
    proxy_set_header Host $host;
    proxy_set_header X-Forwarded-Proto $scheme;
    proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;
}

location /geoserver/ {
    proxy_pass http://127.0.0.1:8080/geoserver/;
    proxy_pass_header Set-Cookie;
    proxy_set_header Host $host;
    proxy_set_header X-Forwarded-Proto $scheme;
    proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;
}

Is there anything I'm missing there?


texts: Going round in Circles..Buffering/ArcObjects/VBA

texts:
I've just started learning VBA and ArcObjects, and at the moment just playing around trying to do different things.
What i am trying to do is Click on a tool, and allow the user to put a point graphic on the map, but also buffer the point(marker) that has been added. I can get the point on the map, but seem to be going round in circles getting the buffer to work. Can someone give me some hints ( i don't want the exact answer just some guidance on where I am going wrong.

Code is here.

Private Sub Buffer_MouseDown(ByVal button As Long, ByVal shift As Long, ByVal x As Long, ByVal y As Long)

Dim pMxDoc As IMxDocument ' Create Object for current document
Set pMxDoc = ThisDocument ' Set object to current Document

Dim pMap As IMap ' Create object for current map
Set pMap = pMxDoc.FocusMap ' Set object to the current active map

Dim pPoint As IPoint 'Create object for point on map
Set pPoint = pMxDoc.CurrentLocation '  Set object to current location

Dim pElement As IElement ' Create object for graphic THAT IS ADDED TO THE MAP!!!!
Set pElement = New MarkerElement ' Set object to a new marker element (point graphic)

pElement.Geometry = pPoint ' Assign the geometry of the element to the Point.

Dim pGraphics As IGraphicsContainer ' Create object to hold the graphic
Set pGraphics = pMap ' Set object to current active map

pGraphics.AddElement pElement, 0 '  Set the graphic that will be added

Dim pActiveView As IActiveView ' Create object for active map display
Set pActiveView = pMxDoc.ActiveView ' Set object to the Active view

'--------------CREATE BUFFER------------------------'

 Dim pSpaRef As ISpatialReference3
 Set pSpaRef = pMap.SpatialReference

 Dim pFCBuffer As IFeatureCursorBuffer2 ' Create Object for buffer
 Set pFCBuffer = New FeatureCursorBuffer ' Set object for new buffer

 pFCBuffer.Dissolve = True
 pFCBuffer.ValueDistance = 100
 Set pFCBuffer.BufferSpatialReference = pSpaRef
 Set pFCBuffer.DataFrameSpatialReference = pSpaRef
 Set pFCBuffer.SourceSpatialReference = pSpaRef
 Set pFCBuffer.TargetSpatialReference = pSpaRef

 Dim pCGLayer As ICompositeGraphicsLayer2
 Set pCGLayer = New CompositeGraphicsLayer

 Dim pGLayer As IGraphicsLayer
 Set pGLayer = pMap.ActiveGraphicsLayer

 pFCBuffer.BufferToGraphics pGLayer

   pActiveView.PartialRefresh esriViewGraphics, pElement, Nothing ' [partial refresh of the graphics element.

  End Sub


texts: How to reproject a latlong J2000 coordinate image into latlong Galactic Coordinate image

texts:
I have a 2x1 ratio lat/long projection of the night sky in J2000 equatorial coordinates (the fixed equatorial coordinate system on Jan 1.5, 2000).

I need to convert this into a 2x1 ratio lat/long image based on the galactic coordinate system (with the center of the image centered on the center of the milky way).

I know that the North Pole of the J2000 coordinate system is at 122.93200023 long, 27.12843056 lat, and 12.860114 East of North.

It's just a rotation and translation of the coordinates, so I would assume GDAL could handle that, but I can't find a function that seems to work.


texts:
I've got a simple problem: I want to count the number of points within a set of polygons.
I have a SQL already but it only gives back the gid of the polygone that actually contains points.
My tables: a polygon layer with 19.000 rows and a point layer with 450 rows.
The following SQL-query
SELECT grid.gid, count(*) AS totale
FROM grid, kioskdhd3
WHERE st_contains(grid.geom,kioskdhd3.geom)
GROUP BY grid.gid;

returns only some 320 polygons that actually contain points. But I want all polygons returned, even thought the number of points is 0.
Of course it has to do with my WHERE-clause. Where do I have to put in my st_contains()?


texts:
I've been trying several times to sample (Sample tool from Spatial Analysis ArcGis 10.5.1) a list of rasters to a location sample (also in raster format .tif).

Here's the code I've been using, and apparently the problem raised up when Python execute the line of the tool - it creates well all the files in the specified folder:

import arcpy, os, sys, string
from arcpy import env
from arcpy.sa import *

arcpy.env.workspace = r'C:\Users\Victoria\Desktop\***\***'
inputFolder = arcpy.env.workspace
outputFolder = os.path.join(inputFolder, "sample")

raster = arcpy.ListRasters('*br*')
colony = arcpy.ListRasters('*col_1*')

for r in raster:
    outName = os.path.splitext(r)[0] + "_col1"
    saveDir = os.path.join(outputFolder, outName)
    Sample(r, colony, saveDir, "NEAREST", "Value", "CURRENT_SLICE")
    print '{0} saved in {1}.'.format(outName, saveDir)

And here's the error that shows up after run the code:

    Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.5\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\Users\Victoria\Desktop\***\***\prueba1.py", line 15, in <module>
    proc = Sample(r, colony, saveDir, "NEAREST", "Value", "CURRENT_SLICE")
  File "C:\Program Files\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 1617, in Sample
    process_as_multidimensional)
  File "C:\Program Files\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 1609, in Wrapper
    process_as_multidimensional)
  File "C:\Program Files\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool

I'm really new in this.


texts: Getting active PostgreSQL connection's username as variable in QGIS?

texts:
Is it possible to get login username from a current PostgreSQL layer connection in QGIS?

My goal is to automatically populate new objects in a layer with an editor login username (not @user_account_name or @user_full_name, because they are not related to a postgres connection).


texts: ArcGIS JS API Popup does not show properly

texts:
This is regarding arcgis js api popups.
when viewed in web the popups appear in good way, but when the browser is changed to mobile, the popups submerge below the screen.
I would like to know how to set popups to right place in mobile also.

I used the same code present in below link
ArcGIS JS API

Rgds,
Amit


texts: Accessing Attribute Table (in composer) with PyQgis

texts:
I have a QGIS (3.44) project containing a composer. In this composer i inserted an attribute table (which i've given the item_id 'tab').

I would like to be able to access this attribute table properties using PyQgis.
I'm using the following code so far :

projectInstance = QgsProject.instance()
projectLayoutManager = projectInstance.layoutManager()
at=projectLayoutManager.layouts()[0].itemById('tab')

The last line returns a QgsLayoutFrame which is a bit misleading for me ...

How can i 'browse' the content of this QgsLayoutFrame to reach my attribute table ?


texts:
I am following the accepted answer for the question: Batch conversion of netCDF to Raster in ArcGIS using Python, obviously, to convert netcdfs. The code taken from the above answer (pasted below) works fine on a single NetCDF. However, if I run the same function on many .nc files by making inNetCDF variable a parameter, ArcMap breaks down. It crashes with ArcGIS 10.1 and hangs with ArcGIS 10.2.2 (on two different machines) after generating about 500 or 600 images.

I suspected that the issue is that ArcMap can't handle this many layers (more than 1000 considering the extracted nc layers). So I added:

      arcpy.env.addOutputsToMap = False

to forbid adding output to the content index, but the crashes still remain.

Is it possible to run the involved ArcToolbox functions like CopyRaster_management outside of ArcMap from command line python and, if not, how do I avoid the crash problem?

--code--

def extractAllNetCDF():

    variable = "RRt_10m"
    x_dimension = "lon"
    y_dimension = "lat"
    band_dimension = ""
    dimension = "time"
    valueSelectionMethod = "BY_VALUE"

    outLoc = "E:/New Folder/"
    inNetCDF = "E:/netCDFFiles/RRt.nc"

    nc_FP = arcpy.NetCDFFileProperties(inNetCDF)
    nc_Dim = nc_FP.getDimensions()

    for dimension in nc_Dim:
        top = nc_FP.getDimensionSize(dimension)
        for i in range(0, top):
            if dimension == "time":

                dimension_values = nc_FP.getDimensionValue(dimension, i)
                nowFile = str(dimension_values)

                #THIS IS THE NEW CODE HERE
                dv1 = ["time", dimension_value]
                dimension_values = [dv1]
                #END NEW CODE

                arcpy.MakeNetCDFRasterLayer_md(inNetCDF, variable, x_dimension, y_dimension, nowFile, band_dimension, dimension_values, valueSelectionMethod)
                arcpy.CopyRaster_management(nowFile, outLoc + nowFile + ".img", "", "", "", "NONE", "NONE", "")
                print dimension_values, i

Here's the outcome using python from the command line.
The script works well except for needing a few imports like arcpy and sys. There was still a crash, but the number of processed images increased 10 fold to about 7000. This time, memory seems to be at issue, with 1.8G consumed before the crash. I am trying the 64 bit background processing as suggested.


texts: How to get around "name 'row' is not defined" error when row in question doesn't exist?

texts:
I am getting a "name 'row' is not defined" error when I run a particular piece of code (see below). This only happens when my input layer has one entry instead of two. So when the attribute query is run to find a second row there is no row to find and hence the error appears.

The aim of my overall script mean this will happen from time to time, what I am looking for is a way to accommodate this error and either avoid this section of script when an error is likely to occur or continue with the script even if this error pops up. I have tried to incorporate a 'try' and 'except' statement into the loop but have been unsuccessful. Would this work? if so how? If anyone has an idea on how to get around this problem please let me know.

import arcpy
from arcpy import env

env.workspace = "C:\MyArcGIS\EDRN_Shp"

arcpy.env.overwriteOutput = True

join110 = "C:\\MyArcGIS\\EDRN_Shp\\join110.shp"
joinField = "LCC_DRN_ID"
query1 ='"FID" = 1'

arcpy.MakeFeatureLayer_management(join110, "feature_layer12")
arcpy.SelectLayerByAttribute_management ("feature_layer12", "NEW_SELECTION", query1)

rows = arcpy.SearchCursor("feature_layer12")
for row in rows:
    if row.getValue("edrn_in") >= "E":
        arcpy.JoinField_management("feature_layer10", joinField, "feature_layer12", joinField, "edrn_in")
    if row.getValue("edrn_out") >= "E":
        arcpy.JoinField_management("feature_layer10", joinField, "feature_layer12", joinField, "edrn_out")

del row, rows


texts: How do I trick QGIS's map composer into chopping off the 000's in UTM coordinates?

texts:
In QGIS 2.8 (ubuntu), you can create a map with ticks and labels. I want the first two digits to be smaller sized, the next 3 to to be normal sized, and hide the last 3 digits. QGIS does not offer this feature directly.

I tried to hide the last 3 digits by creating a custom projection. I just changed +units=m to +units=km. This seemed to work in the main screen, but the CRS behaves like the original CRS in map composer. I don't know if this is user error or a bug in QGIS. Maybe I just have to dive further into custom CRS and change 1.0s to 0.001s (or the other way) until it works?

Furthermore, I would like to get just the first 2 digits. If I run cs2cs -lu, a kilometer is the longest metric distance. Can I add a unit for 1000 kilometers? Or do I need to dig further into the sub-definitions?

Finally, is it possible to define a coordinate system that just displays the middle digits? Instead of dividing by 1000, I want to get the remainder instead.


texts: Stretching only the Y dimension in polygons

texts:
I have a shapefile created in another program, when imported into ArcMap I found that the y axis has been squashed by a factor of 5.

What would be the easiest way of stretching it out, by a multiple of 5 without affecting the X axis?

These are complex polygons, and sorting this out in the source program is not possible.


texts: Removing Dashed Borders of CartoDB Polygons

texts:
I've created a CartoDB map with a bunch of different multipolygons and it is displayed with a CartoCSS color ramp I created.  However I've noticed that I cannot rid my map of the dashed lines that separate the polygons. I initially thought this was due to my polygons being mulitpart instead of singlepart, but even working with a singlepart polygon, I still get those dashed lines.

My CartoCSS has no mention of line-color, line-width, or line-opacity to it.  I've also tried to specifically set the line-opacity and line-width to 0 in hopes that this would take care of it. Here is my CartoCSS:

#test{
  [fcst<=68]{polygon-fill:#DFF22C;}
  [fcst<=66]{polygon-fill:#BFE526;}
  [fcst<=64]{polygon-fill:#9FD81F;}
  [fcst<=62]{polygon-fill:#7FCC19;}
  [fcst<=60]{polygon-fill:#5FBF13;}
  [fcst<=58]{polygon-fill:#3FB20C;}
  [fcst<=56]{polygon-fill:#1FA506;}
  [fcst<=54]{polygon-fill:#009900;}
       ....etc.....
  [fcst<=6]{polygon-fill:#660066;}
}

Here's a screenshot of what I'm seeing

Here's a URL to my visualization: http://cdb.io/1Ig5cGq

It looks like what I'm seeing is actually slivers of the basemap behind the polygons. @andytilia's comment helped me find this. This is interesting because the polygons match up completely.  Even when zooming in all the way the "gap" only appears to be a very faint dotted line.

Here's what I see on a dark map background:

Zoomed in all the way:

I'm wondering now if this is an bug in CartoDB or it's rendering engine Mapnik.  Or perhaps there is something I need to do to my geometry (ST_Scale?) that I could use to "grow" my geometry by like a half meter and remove those lines?


texts: When to rescale reflectance data to calculate vegetation indices

texts:
I'm attempting to calculate a number of vegetation indices using hyperspectral AVIRIS data in ENVI. Some of these indices are normalized, such as NDVI, whereas others are not, such as the Transformed Chlorophyll Absorption Ratio (TCAR).

For the normalized indices such as NDVI, not doing anything to the bands produces an output raster within the valid range. However, in trying to calculate something like TCAR, not doing anything to the bands produces values well outside the valid range.

I have not been able to find anything definitive on when to rescale reflectance data when calculating vegetation indices. Am I correct that all normalized indices can take the unaltered reflectance data, while all other indices must have the reflectance data rescaled by multiplying the individual bands by 0.0001? In particular, is it true that indices that have independent parameters, such as the soil-adjusted vegetation index (SAVI), must also have reflectance data rescaled by 0.0001?


texts: Changing global template features to reflect in other maps

texts:
I'm working on a project that involves 20+ map figures.

We currently are working with a .mxd in which we are using as a template for these figures.

It's most likely that our client will want some global edits made that will need to be reflected in the template (i.e. street layer color changed, or removed element).

Is there an easier way (other than to make 20+ individual global changes to each figure) to have template changes reflected in all of our figures?


texts:
I have a polyline as GeoJSON layer in Leaflet.

I need a JavaScript event when the user click on the portion of the line.

I have to highlight the selected line segment. Is it possible this?


texts: arcpy.env.qualifiedFieldNames = False not working when using arcpy.AddJoin_management

texts:
I would like to add a join on a layer while keeping the original fields names. I thought setting the qualifiedFieldNames environment parameter to False would do the job, as the documentation says, but it has no effect.

import arcpy, os
arcpy.env.qualifiedFieldNames = False

input_fc = os.path.join("C:\\", "temp", "GEODEV.sde", "dataset", "my_feature_class")
join_table = arcpy.TableToTable_conversion(os.path.join("C:\\", "temp", "to_join.csv"), "in_memory", "join_table")
layer = arcpy.MakeFeatureLayer_management(input_fc, "in_memory\\layer")
arcpy.AddJoin_management(layer, "MY_IN_FIELD", join_table, "MY_JOIN_FIELD")

print([f.name for f in arcpy.ListFields(layer)])

>>>'["GEODEV.dataset.field1", "GEODEV.dataset.field2", "GEODEV.dataset.field3", "join_table.fieldA", "join_table.fieldB", "join_table.fieldC"]'

What am I doing wrong?

T.h.a.n.k you a.l.l for your h.e.l.p!


texts: Selecting new records using values inside previous selected rows in QGIS?

texts:
I am trying to create a selection. First, I select records which have Null value in "Status" column ("Status" is NULL"). Then it starts to be a bit tricky. I want to select all the records which have the same number as selected "Status" records. Is it possible to do it in QGIS?

First, I select null values using "Status" is NULL expression.

Then I need to add something to take ID values and use it to create a new selection.

The final selection should look like that:

I am using QGIS 2.16.


texts: How in-sync is the MapBox flavour of Leaflet with the FOSS version?

texts:
How in-sync is the MapBox flavour of Leaflet with the open source version?

I'm concerned that if I got with the MapBox sku, I'll be locking myself into a non-standard version, so things like ... are worrying me

if there is a new feature in Leaflet, how soon does it appear in MapBox ... assuming it does?
do plugin's work with both flavours?


texts: Video Game Controller as a mouse to use in GIS

texts:
I work using ArcFM mapping electric utilities and am rather new to GIS, my question is would it be possible to hook up a video game controller like an Xbox controller to my computer and use it as a mouse to use in GIS? Has someone tried this before? Has it helped in using GIS more efficiently?

I love using my keyboard to automate certain commands that are repetitive and was thinking if maybe something like a video game controller I could automate certain commands to certain buttons that would help me do work faster.

The purpose of this post is to know if other GIS users have already used a different controller other than a mouse and to know if it was better or more comfortable to use GIS with a different controller


texts:
I was wondering if there is a list of steps I can take to screen my variables
for assumption violations prior to using them in my GWR. I have watched this video
tutorial :

https://www.youtube.com/watch?v=plfCMZhROeQ
and read this post:

Regression Model and Spatial Autocorrelation

However, I am still a little bit confused. I understand from logic and background information the types of explanatory variables I expect to explain my dependent variable (these account for about 4-7 from the >20 exp. variables I have. What I don't understand is how to convince myself I should not use certain ones.

The video tutorial suggesting using Scatterplots and OLS (ex VIF value, non-stationarity etc) prior to using GWR. It was my understanding from readings on this forum that I should not use OLS to screen variables for GWR. Thus, I am a little confused.


texts: ArcGIS Online - Display TMS with inverted Y

texts:
In ArcGIS Online, I'm trying to display a Tile Layer (Add > Add Layer from Web > A Tile Layer) that has the Y (row) axis inverted. I tried placing a minus (-) but it didn't worked, the layer disappeared.
Here are the syntaxes I tried :

https://URL/{level}/{col}/{row}.png   >  Working but row (y) is inverted.
https://URL/{level}/{col}/-{row}.png  >  Not working (nothing displayed)
https://URL/{level}/{col}/{-row}.png  >  Not working (nothing displayed)
https://URL/{level}/{col}/%2D{row}.png  >  Not working (nothing displayed)
https://URL/{level}/{col}/{%2Drow}.png  >  Not working (nothing displayed)

(%2D is the URL encoding for - according to https://www.w3schools.com/tags/ref_urlencode.ASP)
Any other idea of syntax or other way to do it ?


texts: Finding facilities on current direction for Network Analyst result

texts:
We use ArcGIS Network Analyst to solve routing problems.

We received a request about filtering facilities (point feature class) along the route result. Also the direction should match. The facilities on the other side of the road should not be displayed.

How can we do it with or without post processing?


texts: Grass7 - r.regression.multi and new raster

texts:
I am on the learning curve with Grass7. I managed to make r.regression.multi work, but there is one thing that would be of high interest for me: it is possible to create a new raster showing the results (residuals) of the analysis in a new raster?
This would enable me to nicely locate the areas with high residuals.


texts:
I need to apply a .qml style to a layer on my map as part of a Python plugin. I gather that this can be done via QgsMapLayer.loadNamedStyle(). The problem I have is that my .qml file is stored within my plugin directory and I can't figure out how to reference the path at runtime.
The plugin is in my user plugin directory (/home/[user]/.qgis2/python/plugins/[plugindirectory]) not the main QGIS plugin directory, so I can't get a reference via qgsApplication.pluginPath().
What's the best way to accomplish this?


texts: Merging Coordinate (long, lat) to a region(polygon)

texts:
At present I have sports grant data, coded by longitude and latitude. This data is broken down by amount received, as well as the intended sport.

I also have various other data, broken down by electoral division (ED).

I would like to do the following. I would like to merge these individual longitude and latitude points into their respective electoral division (polygons). Thus I would know the particular ED grant recipient A was in, how much they received, what sport, along with all the relevant information about that particular ED (population, employment etc.).

Edit - As the image shows, I am running into some issues

Update 25/2/16 Images follow comments below


texts: shp2pgsql-gui not work after Ubuntu update

texts:
I'm using Ubuntu 14.04 with postgres-9.5 and postgis 2.2.2.

Previously I performed shp2pgsql-gui to access the interface, but now, the terminal shows me this message:

the program shp2pgsql-gui not installed in the moment. you can install
typing sudo-apt-get install postgis

The result of:

select * from postGIS_full_version()

is:

POSTGIS="2.2.2 r14797" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER

Does anyone know how I can solve this problem?


texts: Why do some Customised Symbols look good in map but disappear from Legend on Export?

texts:
I am using ArcGIS Desktop 10.0

I created some customised symbols for my map which ArcGIS Desktop did not have. To do this I used BMP files and imported the symbols.

They look good on the map but do not appear well in the legends when I export the map into any format e.g JPEG.

They either do not appear or appear very light.

What can cause this?


texts: Multi/hyper-spectral file format that stores band information and meta-data

texts:
I usually use ERDAS LAN file format drivers to create a multispectral file that contains multiple rasters of different spectral bands. However, I can't find a way to store the information regarding the spectral band acquired by each raster, and when I open the file using, for example, QGIS I can see only the number  of the band (e.g. Band 1, Band 2 etc. etc.) and so it is hard for me to understand which spectral range I am currently using for false colour composting.

I wonder if there are (possibly open-source) file formats allowing storing multispectral raster and metadata for each spectral band collected.


texts:
I have a list of csv files "name.csv" that I want to use to make various analyses in QGIS. I would like to automatize the following steps for each file :

create new project name_main.qgis;
import text layer;
import polygon layer;
[... some analyses...]
save project to a new directory

I know it is possible to do batch analyses with the Atlas plugin, e.g. to export various maps. I also looked at the processing plugin but it seems to work only within one project. But in this case the different files correspond to different projects, not maps.

Each dataset corresponds to 1 site within 1 region.
I would like to be able to visualize each map, if possible. Hence various projects. But to be more precise, I need to calculate mean slope from polygons (so, raster interpolation, terrain analysis and zonal statistics). I also need to reproject all the maps (coordinates x,y,z in meters) with GPS coordinates I have for 4 points each time.

Maybe I should first reproject from gps points and then merge all the sites in one map to do batch analyses?

I got the first steps to work but I have an issue for importing shapefiles: I need to first create shapefiles from text files, then to re-import the shapefile.

Here is what I have :

for cname in glob.glob('*.txt'):
      uri2 = "file:///" + csv_path + cname + "?delimiter=%s&wktField=%s" % (";", "wkt")
      name = cname.replace('.txt', '')
      cont = QgsVectorLayer(uri2, name, "delimitedtext")
      shlayer =  QgsVectorFileWriter.writeAsVectorFormat(cont,csv_path + name +'.shp',"utf-8",None,"ESRI Shapefile")
# Add to QGIS canvas
iface.addVectorLayer(shlayer, name, "ogr")

The shlayer is not created, and obviously I get an error when I try to import it.
Is there a problem with how I import the WKT layer ?


texts: Crop/clip raster with shapefile that do not overlap

texts:
I have shapefile that intersects with a raster but do not fully overlap :

I would like to clip the raster with this shapefile but every time I try I  get this error message:

ValueError: Input shapes do not overlap raster.

I have tried two methods to clip the raster:
once with rasterio:
#get coordinates of shapefile
def getFeatures(gdf):
      """Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
       import json
       return [json.loads(gdf.to_json())['features'][0]['geometry']]

coords = getFeatures(plots)

out_img, out_transform = mask(img,coords, crop=True)

and once with earthpy:
crop=es.crop_image(img,plots)

both give me the same errors.
This is how I read the shape and the image:
#load shapefile
poly = geopandas.read_file(r'very\nice\path\and\polygon.shp')

#determine projection for shape and change it into the image projection
poly.crs = {'init' :'epsg:4326'}
poly=poly.to_crs({'init': 'epsg:32721'})

#Load image
img=rasterio.open(r'what\a\wonderful\path\image.tif')

after this I used the methods that demonstrated.
How can I clip the raster with the polygon even if they are not the overlap 100%?
Edit:
I have tried this:
out_image, out_meta = es.crop_image(img, plots.geometry)

and also the same with rasterio as was suggested here but for some reason it clips only for the extent of the shape but not to the actual shape:


texts: GDAL: Change no-data value of GeoTIFF without having to recreate the whole GeoTIFF

texts:
I've got a problem with some big GeoTIFFs, that they got no-data value -128, being unsigned. And because of that I can't upload them to GEE, as I am getting an error:
Error: No data value -128.0 cannot be applied to band #0 of type Short<0, 255>.

Is there an automated way to change these no-data values to something else (basically these files are having max 10, so 255 would be fine I guess; or even 0) without having to recreate the whole datasets (that takes ages, as it's the whole Europe in 10m).
Could be in Python or GDAL tools.
Example gdalinfo of a file (with lower resolution, but having the same problem):
Size is 11282, 6117
Coordinate System is:
PROJCRS["ETRS_1989_LAEA",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Lambert Azimuthal Equal Area",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",52,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",10,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",4321000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",3210000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3035]]
Data axis to CRS axis mapping: 1,2
Origin = (-1217000.000000000931323,6415000.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1217000.000, 6415000.000) ( 82d45'43.92"W, 41d20' 9.92"N)
Lower Left  (-1217000.000,  298000.000) ( 40d 5'20.45"W,  9d20'13.32"N)
Upper Right (10065000.000, 6415000.000) (103d45' 8.42"E, 39d27' 3.98"N)
Lower Right (10065000.000,  298000.000) ( 61d46'28.34"E,  8d 9'45.42"N)
Center      ( 4424000.000, 3356500.000) ( 11d32'42.50"E, 53d18'24.09"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  Min=0.000 Max=10.000
  Minimum=0.000, Maximum=10.000, Mean=0.206, StdDev=1.248
  NoData Value=-128
  Metadata:
    STATISTICS_MAXIMUM=10
    STATISTICS_MEAN=0.20599296116556
    STATISTICS_MINIMUM=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=1.2477841778992
  Image Structure Metadata:
    PIXELTYPE=SIGNEDBYTE


texts: I'm looking for a set of outlines of the river basins in the UK - I don't need a particularly detailed or highly accurate map - just a map saying point X will drain into the Thames or point Y will just go into the sea!  Just the main river names would be fine.  I have the full Environment Agency rivers set and the OS DTMs if I want... but how do I convert this into a catchment boundary?


texts: Finding number of Points within Cluster by different sizes using QGIS?

texts:
I'm trying to figure out how many points are within a pre-determined cluster size. For instance, I can create this map

I can see there are a few clusters and I can create a heatmap to illustrate the density. However, I'm trying to figure out how to set a cluster with a radius of 5, 10, and 20 miles and count how many points are within those cluster.

Do you know how to implement this in QGIS?


texts: Compute centroids and buffers for multiple polygons in Google Earth Engine

texts:
I have a dataset with many polygons that come from a fusion table ('xxxxx'):

var polygons = ee.FeatureCollection('xxxxx');

I want to compute the centroid and a buffer of Z meters for each of these polygons. How do I do that?


texts: Constructing a "dumb" lat-long to lat-long conversion in pyproj

texts:
I have:

import pyproj
c = pyproj.Proj(init='epsg:4326')
c(180, 180)

I expect:

(180, 180)

I get:

(1e+30, 1e+30)

Why is that?


texts: Auto-Sanitizing illegal PostGIS WGS84 (lat, lon) coordinates

texts:
I have just stumbled upon what seems to be a very curious behaviour by PostGIS.

The XYZ point -117.0, 32.0, 0 in WGS84 (SRID 4326) is legal:

items=> SELECT st_y(st_geomfromtext('POINT Z (-117.0 32.0 0)', 4326));
 st_y
------
   32
(1 row)

Now consider the XYZ point 32.0, -117.0, 0. It should be illegal under WGS84 (SRID 4326), because earth latitudes should be within `[-90, 90], but PostGIS accepts it:

items=> SELECT st_y(st_geomfromtext('POINT Z (32.0 -117.0 0)', 4326));
 st_y
------
 -117
(1 row)

PostGIS will even calculate the distance between the two illegal points, giving a minor notice:

items=> SELECT st_distance_sphere(
        (st_geomfromtext('POINT Z (32.0 -117.0 0)', 4326)),
        (st_geomfromtext('POINT Z (-117.0 32.0 0)', 4326)));
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
CONTEXT:  SQL function "st_distance_sphere" statement 1
 st_distance_sphere
--------------------
   15938523.6399401
(1 row)

This caused a very nasty bug with an external library which sanely refused to accept a -117.0 latitude.

Can I block PostgreSQL from accepting coordinates violating WGS84 into the DB?


texts: How can I reproduce the ArcMap "Select Layer By Attribute" tool as a Web Service Geoprocessing Tool?

texts:
I would like to reproduce the Select By Attributes tool below in ArcMap as a geoprocessing tool to be consumed as a service in a Portal for ArcGIS Web App.

I am not a developer.

How do I get the sql query builder included in the Select Layer By Attribute geoprocessing service?


texts: Expedite Near function in R for 10^9 coordinates?

texts:
Goal
I want to find the geodesic distance between the closest hydrography feature to each of my lat/long coordinates.

I have hydrography data from New York state in a shapefile. I have many points, 10^9 number of coordinates in WGS1984 datum. There's about 20,000 features in my hydrography data which was originally in NAD1983.

I'm trying to find a way to expedite my code. This is what I've tried:

1) dist2Line from the package(geosphere): I keep getting the error: "returning Infno non-missing arguments to min". When I can get it to work (I haven't figured out why sometimes it works and sometimes it doesn't), it takes too long as well.

2) gDistance from package(rgeos): this also expects planar coordinates and takes too long (~43 sec for 300 pts).

3) st_nearest_feature from package (sf): This is fast (1 sec; 300 pts) but assumes the coordinates are planar. I get the index of the nearest hydro feature and I don't really know how to calculate distance easily from there (I suppose I could find a centroid for each hydro feature and calculate distance to that feature)?

4) st_nearest_points from package (sf): this takes too long (~100 sec, 300 pts) and assumes that coordinates are planar. It gives me linestring so I suppose I could calculate the distance between the two coordinates using a basic distance formula but the time & projection is not ideal.

I'm at a lost at how to proceed in a timely manner from here.

Potential Thoughts:

I've converted all of my hydrography features to WGS 1984 to work with it. Should I consider converting both files into planar coordinates? I'm not sure I understand well enough whether converting from unprojected to projected would cause any potential problems.

I just want the fastest way to process all of these points. I've been using fast_extract from package(velox) and (prioritizr) with raster data, and it's been incredibly fast. Something like that would be ideal.

I've tried to expedite this in ArcMap or R, and am open to either method. I've created a second question about ArcMap in the event somebody can think of a better me

Expedite Near function in ArcMap 10.7 for 10^9+ coordinates?


texts:
I am new to GIS world and trying to incorporate OpenLayers with GML or SVG.

I am building a map engine and one of the features should be support of indoor maps. It can be accomplished with a static image, but I got a requirement that the engine has to understand vector format, the main reason is to keep a good quality when you zoom in or zoom out.

  Common use case: the client has an indoor map (floor plan) made in
  AutoCAD, so it can be exported to GIS formats.

I do not understand how to make it work. I should define the projection but I have no clue which number should be assigned if I only want to show my floor. Regarding the static image, I could create a custom projection, but how to deal with GML?

To be honest I can do it with SVG if someone helps me to correctly set up the image size because now it shows the only limited part of the image.

I read SVG data

 viewBox="0 0 1024 768"

and set extent accordingly

var extent = [0, 0, 1024, 768]

besides that, I pointed out the image size

 imageSize: [1024, 768]

it seems to be a trivial task, but not for me.

I prepared demo (wait a few seconds until it loads  SVG file)

UPD:  I managed to make it work with SVG, but anyway it looks blurry when zoom in and image located not in the center


texts: How to find the geographic coordinates of a pixels of known value in a set of images?

texts:
I am trying to know the geographic coordinates of pixels whose pixel value is known from its corresponding GeoTIFF image/array of interest.

clip = [[[1.0 2.3 3.3
          2.4 2.6 2.7
          3.4 3.2 3.9]]]
       [[[3.0 4.3 7.3
          2.7 2.6 5.7
          3.4 4.2 8.9]]]

#interested pixels whose geographic coordinates are needed from their corresponding image. First value from first image and second value from second image...

values = (3.2, 2.6)

Here, the clip stores all the arrays of interest. The shape of each array of interest/subset image is (1,3,3).

I got stuck here. I am unable to find the geographical coordinates of pixels of interest within these arrays of interest, the pixels of interest are given in values (1st value for the first array and so on).

The output should be the coordinates of 3.2 from the first array/image and 2.6 from the second array.

Expected outcome:

# longitude and latitude of points 3.2 from the first image and 2.6 from the second image stored in the "clip".

Coordinates = (64.3, 18.6)
              (97.2, 11.2)

Any help in this regard is much appreciated.


texts: Pull out multiple elevation points from an ISurface

texts:
I have an ISurface object and wish to pull out several elevations with one call; instead of using the GetElevation(IPoint) method.

The ISurface is representing a dem with 30m cell size.

Is there an ISurface method that pulls out several elevations with one call?

Is there another way of pulling elevations from a dem?  A query?