texts: Setting constraint so that new value in "Table A" must already exist in "Table B", else return error?

texts:
I'm trying to set up some kind of constraint so that I can only enter values in "Table A", if that value exists in "Table B".

I want this to act sort of like a selection list, where there are only a few values that exist in "Table B", which then restricts the values that can be entered into "Table A". (Kind of like a drop down list/ selection list)

I'm using pgAdmin 4 to access a PostGIS DB


texts: How does Google Map API authenticate users?

texts:
They put a limit on the amount of map views, geocodes etc that can be done each day but in the documentation I can't find anywhere where they actually list how it is authenticated.

My guess is it has to be one of the following:

Website url
Ip address
Google account

If anyone knows, that would be really useful.


texts: Statewide (and more) Address Point data source

texts:
I am seeking a data source for a statewide address point layer. I do not need resident or business information at these addresses; I just need the validated structured addresses (preferably with standardized elements) and the coordinate.

This is for an address suggestion application on clients without internet access. We could do server side queries to an internet service, but not client side. (So, for example, this is why we cannot use Google for queries.) We will be querying partial addresses to build suggestions, so generally geocode services will not work. Partial addresses will have a zip code.

Ideally though, we would obtain a database of address points, allowing a partioned table query server side against that database.

The main State we want is Missouri, followed by Illinois, but could go as high as a nationwide database. We have interest at this time in countries other than the United States.

Do not worry about budget issues on this. I need to gather information on available sources before working up a budget/RFP.


texts: Stack two numpy array into one image with two bands in order to clip it with anotherr shapefile

texts:
I have two np array,'array 1' and 'array2' , and I have one image with 3 bands.
I want to stack array 1 and array two to be one image with two bands with the same shape as my image and then to clip it with another shapefile that I have.
Array 1 and Array 2 are coming originally from the image so there are enough pixels to be the same shape.
This is how I have tried to do that:
#Load original image and check its' shape:
img=rasterio.open("original.tif")
array=img.read()
array.shape

>>>(3,2199,4041)

Here I did many calculations that in the end created one pandas df when each pixel from Original image is a row withour the shape of original image. I have created two numpy arrays from this final table which I want to construct into one image with to bands in the end:
Array_one=df_matrix['one'].values.reshape(2199, 4041)
Array_two=df_matrix['two'].values.reshape(2199, 4041)

Then I stack the two arrays together:
stack=np.stack((Array_one,Array_two))

stack
>>>array([[[0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        [0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        [0.6, 0.7, 0.7, ..., 0.2, 0.2, 0.2],
        ...,
        [0.5, 0.5, 0.5, ..., 0.2, 0.2, 0.3],
        [0.5, 0.5, 0.5, ..., 0.2, 0.2, 0.3],
        [0.4, 0.4, 0.5, ..., 0.3, 0.2, 0.3]],

       [[0.3, 0.3, 0.3, ..., 0.3, 0.4, 0.3],
        [0.3, 0.2, 0.2, ..., 0.4, 0.4, 0.3],
        [0.2, 0.2, 0.2, ..., 0.3, 0.3, 0.2],
        ...,
        [0.3, 0.5, 0.5, ..., 0.2, 0.2, 0.2],
        [0.4, 0.4, 0.5, ..., 0.2, 0.2, 0.2],
        [0.3, 0.4, 0.4, ..., 0.2, 0.2, 0.2]]])

Now things start to get complicated. I wanted my numpy to be saved as raster so i'll be  able to clip it with rasterio or gdal but I couldn't find a way to do it. also to visualize it failed all the time:

# Use bilinear interpolation (or it is displayed as bicubic by default).
imshow(stack, interpolation="nearest")
show()

TypeError: Invalid shape (2, 2199, 4041) for image data

that also happens if I use show(1).
so then I tried to save it as raster with two band so I can load it again and then clip it , kike this:
 with rasterio.open('test',
                    'w',
                    driver='GTiff',
                    height=stack.shape[0],
                    width=stack.shape[1],
                    count=2,
                    dtype=stack.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(stack, 1)

but it fails all the time with this error:

ValueError: Source shape (1, 2, 2199, 4041) is inconsistent with given
indexes 1

I pretty much sure that there is something i'm missing regard the number of bands.
My endgoal is to create raster which has two bands- array1 and array 2, so I can clip it with my shapefile.
Edit:
I have tried the solution suggested by urban87 , I could display the numpy array :

but when I saved it as raster and load it again like this:
 with rasterio.open('saveThis.tif',
                    'w',
                    driver='GTiff',
                    height=stack.shape[0],
                    width=stack.shape[1],
                    count=2,
                    dtype=stack.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(stack[0], 1)
        dst.write(stack[1], 2)

img=rasterio.open("saveThis.tif")
show(img,1)

I got this raster:


texts: Looping through lists with arcpy.AlterField_management?

texts:
I can't get my code to run correctly and I'm sure have the logic wrong. I have two lists declared as parameters: one is a list of field names that I would like to change from and the other is a list of field names I would like to change to.

I can get the python to run correctly if I first check to ensure the table has all of the old field names and does not have any of the new street names. But, what I would like to do is automate that 'check' - so that the script first checks to see if the old field name exists and, if so, continues with the AlterField task. If the old field names does not exist, it skips that particular change and continues on to the next one.

Here is the (bad) code that I have so far:

import arcpy
arcpy.env.workspace = r'C:\data\MyGDB.gdb\Addresses_WGS84'

fcTable = 'Copy_of_roads_from_sde'
oldfields = ['STR_NAME1', 'SPEED_MPH', 'LZIP', 'RZIP' ]
newfields = ['ST_CONCAT', 'SPEED_IMP', 'ZIP_L', 'ZIP_R']

for field in oldfields:
    if field.name in arcpy.ListFields(fcTable):
        arcpy.AlterField_management(fcTable, oldfields, newfields)

I've tried a few different approaches, stumbling around trying to use break, elif, and else statements in various places, but to no avail. The error I keep getting (when I actually get the syntax right) is:

Failed to execute. Parameters are not valid.
ERROR 000728: Field STR_NAME1 does not exist within table
Failed to execute (AlterField).

I'm sure that I have the overall, big-picture concept wrong. Also, FYI: my code continues on, after altering these field names, to use AddField_managment to add a bunch of other, new fields to my data.


texts: How can I use a variable to specify which field to update with arcpy.UpdateCursor?

texts:
I am writing a script that cycles through several different point shapefiles, finds how many lie within a given area and the writes the outcome to a table in an Access database. I've used an search cursor to look up the point file names, which is then saved to a variable. I'd then like to use this variable to specify which field to update in another table. Here is my code:

selection = 1
while selection < 3:
        rows = arcpy.SearchCursor("C:\\MyDB.mdb\\Source_Class_combinations") #opens search cursor on the info table
        for row in rows:
            if row.Order == selection:
                destfield = row.DestField # extracts the name of the destination field
        del row
        del rows
        print "Defining variables..."
        ptfile = str(destfield)+"_lyr"
        rowtoud = "row."+str(destfield) #create the variable with the field to update

        # Process: Select Layer By Location
        print "select by location"
        arcpy.SelectLayerByLocation_management(ptfile, "INTERSECT", neighfile, "", "NEW_SELECTION")
        count = int(arcpy.GetCount_management(ptfile).getOutput(0))
        print "count = " + str(count)

        #Update the number of points selected in the neighbourhoods table
        print "updating count..."
        rows = arcpy.UpdateCursor("C:\\MyDB.mdb\\neighbourhoods")
        for row in rows:
            if row.StudyID == str(neigh):
                rowtoud = count
                rows.updateRow(row)
        del rows
        del row
        selection +=1

The selection all seems to work fine, but the table is not updated and no error message is returned. Can I use a variable like this?

I'm using ArcGIS 10, python 2.6, and Access 2007

Thanks, Flo


texts: QGIS @parent or import to database(Error 7) not working for layer

texts:
I have a shape file layer that won't load into the SpatiaLite database and it gives Error 7.  It also won't run functions that work in other layers and works fine without the @parent or if I remove the filter.

> aggregate ( 'Mainland_Address_079da3c1_0048_4cf5_ba0d_e89a4d130688',
>              'concatenate',
>              to_string("CIVICNUM"),
>              filter:=intersects( $geometry, geometry(@parent) ) ,
>              concatenator:=' ,')

I can't figure out what is wrong with it.  Unfortunately I can't share it  because it has sensitive information. But I can share the Information and Source Fields.  My project is using EPSG:4326. What is the issue with the layer?

Edit: Error 7 Picture

Edit: Table picture


texts:
For a personal website I would like to plot coordinates of for example historic events etc. on an existing image.

This can be a photograph of an area of which I look up the coordinates of the 4 corners (put those in the database with the url of the image).

Based on the coordinate of the historic event I can select which photographs contain that coordinate and draw the event-coordinate on the selected picture.

This image should be created in php. What is the best solution to do this?

Sorry I'm a newbie in GIS, I only have this idea which I'd like to execute and so I came here.


texts: ArcPad Jscript ThisEvent.Object, autofill form with multiple pages

texts:
I have a script that I'm trying to use to autofill the controls of an ArcPad form with controls on multiple pages.  The 'autofill()' function is referenced in the form's onload event and I would like to avoid opening each page individually.  The following code works well for the first page but the second page does not autofill.  Am I referencing the second page incorrectly?  I have a third page to add as well.  The myentries[] is an array that references specific values for each control (which are comboboxes).

Code:

function autofill(){
var pages, datepage, controls1, notespage, controls2;

datepage=ThisEvent.Object.Pages("date");
controls1=datepage.Controls;
controls1.Item("cboemp").ListIndex=myentries[0];
controls1.Item("cboaction").ListIndex=myentries[1];

notespage=ThisEvent.Object.Pages("notes");
controls2=notespage.Controls;
controls2.Item("cbolarvae").ListIndex=myentries[2];
controls2.Item("cbopupae").ListIndex=myentries[3];
controls2.Item("cboadults").ListIndex=myentries[4];
}


texts: How to switch to a GDB created in the middle of a script

texts:
EDIT:
Here's my full script and hopefully a better way of explaining what I need done:

#Import

import arcpy, sys, os

# Set environment settings

arcpy.env.workspace = sys.argv[1]
inFiles = sys.argv[2]
outWorkspace = sys.argv[3]
projectGDB = sys.argv[4]+".gdb"
coordinate = sys.argv[5]
clipBoundary = sys.argv[6]
scratchGDB = sys.argv[7]+".gdb"

# CREATE A GEODATABASE

arcpy.CreateFileGDB_management(outWorkspace, projectGDB) # Here I create the database
# for the projected files to go into
arcpy.CreateFileGDB_management(outWorkspace, scratchGDB)

# CHANGE THE PROJECTION

arcpy.BatchProject_management(inFiles, outWorkspace+"\\"+gdbName, coordinate)

# CLIP FEATURES TO STUDY AREA BUFFER

fcs = arcpy.ListFeatureClasses() # Here I need it to somehow read from the GDB with
# my projected files in it that was created above..however it reads from the workspace
# and because the workspace is set before running the script I can not choose the GDB
# because it is not created yet

for fc in fcs:
    arcpy.Clip_analysis(fc, clipBoundary, outWorkspace+"\\"+scratchGDB+"\\"+fc+"__clipped")

I am wondering if there is a way to change the workspace as you move through the python script without having to do something like wrap tkinter around it.
For example:
I have set up my inital workspace through the parameters (sys.argv) so that I can create a gdb to store my files - this is a general "Project" folder as my workspace
Next I have a batch project in which the files that the user wants to be projected are chosen, the output of these projected files is to the geodatabse that was created
NEXT (and here's where my issue comes in) I need to clip my projected features to a boundary of the study area, so I have a batch clip written out:

fcs = arcpy.ListFeatureClasses()
for fc in fcs:
     arcpy.Clip_analysis(fc, clipBoundary, outWorkspace+"\\"+gdbName+"\\"+fc+"__clipped")

However to select the fc from fcs I need it to be looking into the projected files GDB to read all the feature classes...however if I just run the script from the beginning to this point the workspace is still sitting in my main project folder and therefore there are no feature classes that will be read into the script........


texts: Creating a script to track down the source image of a pixel within an imagecollection in GGE

texts:
I have created a script of sentinel-2 images from 2017 oct to 2018 april with 20% cloud-cover at 15-days time interval. I have one imagecollection every 15 days, i also calculated NDVI and added that as a band for the first 15-days collection.

Now, every collection has a number of images with specific NDVI values. What i want to do is to compare pixel by pixel NDVI values every image within the collection and choose the pixel with highest NDVI value and print down which image among all the images is giving the highest NDVI and choose the reflectance values corresponding to that pixel and that image within that imagecollection, this way create a 15-day maximum reflectance composite derived from the maximum NDVI value-showing image inside an imagecollection. So, my aim is to get maximum reflectances of images within imagecollection for each pixel. I want to do that for every 15-days collection.

I am very much new to coding but i am in my thesis now and stuck at this point. Can anybody tell me how do i do this in GGE? I am mentioning the first 15-days code here-

var tanguarhaor = ee.Geometry.Polygon([
  [[90.9952, 25.16], [91.2616,25.1824],[91.4841,25.1053],[91.55, 24.9647], [91.4594, 24.8065], [90.8772, 24.814],[90.781,24.9635],[90.8043,25.1103],[90.9952,25.16]]
]);

Map.centerObject(tanguarhaor);

Map.addLayer(tanguarhaor, {color: 'FF0000'}, 'geodesic polygon');

function maskS2clouds(image) {
   var qa = image.select('QA60');
    var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
   var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
              return image.updateMask(mask).divide(10000)}

var start = ee.Date('2017-10-01');

var finish =ee.Date('2017-10-15');

var TahirpurOCT17first = ee.ImageCollection('COPERNICUS/S2')

.filterBounds(tanguarhaor)

.filterDate(start,finish)

 .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',10))

 .map(maskS2clouds);

var IMAGE1 = TahirpurOCT17second.map(function(image) { return image.clip(tanguarhaor); });

var addNDVI = function(image){var ndvi=image.normalizedDifference(['B8','B4']).rename('NDVI');return image.addBands(ndvi);};

var ndvi1 = IMAGE1.map(addNDVI);

var ndvi2 = ndvi1.sort('system:time_start', false);

var collection = ndvi2.select(['NDVI']);

print(collection);

I am stuck at the point how do i know from which image the maximum NDVI value is coming from within the imagecollection 'collection' in the script for each pixel and get the reflectance values corresponding to that pixel of that image within that 15-day imagecollection? Is it possible to get the result as a map somehow?


texts: How do I calculate curvature using GIS co-ordinates

texts: I have GIS data for which I would like to calculate where curvature is present. Anyone has done this before, please help.


texts: How to convert 32-bit raster to double precision?

texts:
I have a raster that is currently 32-bit float, but need to convert it to a data type of Double Precision in order for it hold larger values within each cell (I'm doing some calculations). How can this be done in either ArcGIS Desktop or ArcGIS Pro?

I have tried using Copy Raster, but only the 64-bit option for data type is available and this doesn't have the output desired.


texts:
We have a serious problem with QGIS 2.18. We do design for telecom and we have recently moved from MapInfo to QGIS. We receive MapInfo TAB files from Client and Open them in QGIS. After several days of design, all of sudden QGIS Acts weirdly and it moves the features to random places.
this picture shows where the feature suppose to be(it was there till yesterday)

this is the new place some of the features moved

As the second picture shows all the features have random(invalid coordinates)

Is it something we are doing wrong in QGIS?


texts: Creating boundary polygon shapefile from set of LAS files using ArcGIS Desktop?

texts:
I have a bunch of LAS file sets that I need to process. The current task at hand is to generate a "footprint" shapefile that contains polygons for each LAS tile. I am looking for a graceful way to generate this polygon shapefile for all the LAS files in a directory. I am using ArcGIS 10 with just about every extension you can have (academic site license).

At the moment, my workaround is to use LAStools lasboundary.exe and then merge all of the individual shapefiles into single shapefile. However, this takes a while to execute and requires the additional merge step at the end.


texts: Geoprocessing Service on ArcGIS Server to Export File

texts:
ArcGIS Server /
Geoprocessing Question

I have a geoprocessing service that will perform a selection on a layer based on a provided expression.  I would like to output the result as a csv file (or other type of file) and have it output to the browser for the user to download.  Currently, i can save the output to a scratch directory on the server but I want to return it back to the user.  This is similar to the clip/zip/ship process.  How do I run a process and return a file back to the client?

Does anybody have a detailed workflow or example to point me to?


texts: Catching data frame caption changed event

texts: Is there a way in ArcObjects to catch the event that the data frame's caption has been changed or a layer's caption has been changed?


texts: QGIS SAGA Channel network error: Input layers do not have the same grid extent

texts:
I am trying to follow along with this exercise using my own DEM.

I have done the first step, creating the catchment area, with no trouble. When I try to run the Channel network algorithm, however, I get the error message “Input layers do not have the same grid extent.”

I don’t see anywhere to specify a change in extent or cell size, so I’m puzzled about why the extents are different at all. And then if they are different, how can I get them to be the same? I’ve tried using the raster layer clipping to make new layers from the DEM and the catchment area to the same extent, but this gave me the same error.

I did try running the Channel network and drainage basins algorithm, and this worked but the channels seemed to be incomplete − it looked like only the top two order channels were there, with nothing below these.

I’m using QGIS 2.6.0 on Windows 8.1.

I’ve been using TauDEM to do this work since I haven’t resolved the issue described above.


texts: External commands not working in Geocortex Essentials 4.4 Workflow?

texts:
These following external commands are not working at all when i try to find a feature set thru SQL non-query and then want them to zoom or highlight in Geocortex Essentials Workflow. (Run External command utility)

HighlightESRIFeatureSet
ZoomToESRIFeatureSet
ZoomToFeatures

The Geocortex SDK API Reference for HTML 5 viewer 2.5 displays the above commands as part of it.

But nothing happens over map. Any idea?


texts: Why is the imported TIFF file distorted with noise in Quantum-GIS?

texts:
I've a TIFF-file. Loading the raster image into QGIS does not display correctly. (Using the Add raster layer button.) It seems to be distorted with noise. See image below:

The source of the TIFF file is a 300 dpi scanned document. Other scanned documents display correctly.

Projection doesn't matter as the image is not georeferenced yet.

I'm using qgis-1.7.0 in a sabayon-linux environment.

What's the source of this problem? Did anyone else encounter this before? Thanks for any hints or experience on this topic.


texts: How to assign a new name for values selected from two columns?

texts:
I'd like to assign a new name (either HIGH, MEDIUM or LOW) to values selected from different columns in ArcMap 10.3

I selected the two values from the two columns ("MAIN_ROCK" and "SUB_ROCKS") using SELECT BY ATTRIBUTES as follows:

"MAIN_ROCK"='gravel" AND "SUB_ROCKS" in ('mudstone tephra loess','sand mud','sand mud peat','sandstone mudstone tephra silt lignite','sandstone siltstone mudstone limestone silt lignite')

I want these highlighted main and sub rocks to be called "HIGH" in new column. So, I created another text field and named it CATEGORY.

I used field calculator while these main and sub-rocks are highlighted and wrote 'HIGH' in the field calculator of CATEGORY column but I got the message below

Any idea how can I do this?


texts: FrameCloud Popups from json file in OpenLayers 2.13

texts:
So my problem is that i get only the last popup message from my json file.
My marker builder works nice.
Any thoughts will be helpful guys. Here is my code:

      function airports() {
                pointsmarks = new OpenLayers.Layer.Markers("Airports");
                $.getJSON('../json/airports.json', function (data) {
                    $.each(data, function () {

                        px = this.FIELD7;
                        py = this.FIELD8;
                        name=this.FIELD2;
                        size = new OpenLayers.Size(30, 30);
                        offset = new OpenLayers.Pixel(-(size.w / 2), -(size.h / 2));
                        airicon = new OpenLayers.Icon('../icons/airportsvg.png', size, offset);
                        map.addLayer(pointsmarks);
                        pointsmarks.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(py, px), airicon));

                       pointsmarks.events.register('mouseover', pointsmarks, function(evt) {
                        popup = new OpenLayers.Popup.FramedCloud("Popup",
                        (new OpenLayers.LonLat(py, px)), null, '<div>'+name+'</div>', null, false);
                        map.addPopup(popup);
                         pointsmarks.events.register('mouseout', pointsmarks,
                          setTimeout( function(evt) { popup.destroy(); }, 4000));
                        });
                    });
                });
            }


texts: Division of a raster into multiple rasters by its value

texts:
After the supervised classification of the original raster I want to divide/split the classified raster into multiple rasters by its classes. My aim is to use the single classes as polygons and calculate their areas.
Is there any tool I can use in QGIS?


texts:
I have a shapefile with polygons and want to add an attribute field with the WKT in QGIS 3.4.

The function geom_to_wkt in the Field Calculator isn't working and I don't know how to do it. Any ideas?


texts: Parameters for ED50 -> ETRS89/WGS84 seven-parameter transformation which are good for Denmark?

texts:
I am looking for parameters for a seven-parameter transformation (Rx, Ry, Rz, Tx, Ty, Tz, s) from ED50 UTM (zone 32U) to ETRS89 and/or WGS84.

The best match I could find is here: http://geodesy.spacecenter.dk/~geod/REF/3drefdk.htm

Unfortunately, it seems like these parameters produce coordinates which are off by some 50 meters or so.

Ideally, the parameters should lead to a transformation within a couple of meters of what the online calculator at http://valdemar.kms.dk/trf/ produces.


texts: OpenLayers3 add map to bottom of a layer?

texts:
According to OpenLayers 3 documentation, the function

  addLayer(layer)

Adds the given layer to the top of this map.

How can I add a layer to the bottom of a map?

P.S. I need this to add a background layer similar to a base layer in OpenLayers 2.


texts: Is ArcGIS API for JavaScript compatible with Internet Explorer 7 and 8?

texts:
I've developed a web map using ArcGIS API for JavaScrip version 3.8 and its not compatible with IE 7 or 8. I've advised the people who have reported this to update their version of IE but, if it doesn't involve rewriting the website, I'd like to offer something that can work with these versions.

Here's the website: http://il.water.usgs.gov/data/mchenry/map.php

As you can see, in IE 7 or 8 the map doesn't even load which throws javascript errors.

EDIT

I have added some script that changes the markup and css to handle HTML in IE versions lower than 9 found here: https://stackoverflow.com/questions/5449502/html5-and-css3-for-ie7-and-ie8

Since doing this, I've noticed that the map loads on zoom in/out but none of my points load.


texts: How can I secure my WMS layer URL on GeoServer as whenever a WMS layer is loaded from GeoServer than the layer URL can be easily seen in the browser console?


texts: Filtering clustered markers on Leaflet.js using slider from nouislider?

texts:
I have clustered markers on the map and I am using a library called nouislider for my slider. I am using this to show certain markers based on the range of the distance specified.

I got this working for the markers but no luck on clustered markers. It is removing all my cluster markers rather than filter it. I'm not sure how to approach this.

Here is a link to my code. I can't get the slider to show up on Jfiddle and also got an error on there -"L.geoJson is not a function"

Jfiddle

this is my filter which is working for markers but not for clustered markers, I think the problem is I am not filtering the clustered markers and I am not sure how to do it.

        var slider = document.getElementById("slider");

        noUiSlider
          .create(slider, {
            start: [min + 1, max - 1],
            tooltips: true,
            connect: true,
            range: {
              min: min,
              max: max
            }
          })
          .on("slide", function(e) {
            console.log("e", e);
            surfSpotsGeoJSON.eachLayer(function(layer) {
              console.log("layer", layer);
              if (
                parseFloat(layer.feature.properties.distance) >=
                  parseFloat(e[0]) &&
                parseFloat(layer.feature.properties.distance) <=
                  parseFloat(e[1])
              ) {
                layer.addTo(map);
              } else {
                map.removeLayer(layer);
              }
            });
          });
      })
      .catch(err => console.log(err.message));
  }


texts: How do I tile raster input to speed up display in R Shiny programs?

texts:
Disclaimer: I have very little GIS experience

I've been working on a project in R that involves displaying large raster maps in an R Shiny interface. I am using the leaflet package for display. The displaying process (grabbing the raster from the stack, removing the old raster and adding the new one to the map) takes 4 or 5 seconds, and I'd like to bring that down as much as possible. The data is quite large, but it comes to me in small pieces, maybe keeping it in small pieces could be helpful.

I've already taken a few steps that speed up the display:

I aggregate the data ahead of time so that I'm not aggregating when I display it
I project the rasters for leaflet ahead of time so leaflet doesn't need to reproject when it displays them
I set my <= 0 values to NA ahead of time (my values should be strictly positive, NA values should be transparent on the map)

This means there's no actual processing when the map is displayed, it's just the displaying itself that takes time.

People often ask why I don't "tile the data like google maps" and truthfully that's because I don't quite know what that means and couldn't figure it out online. To the best of my knowledge it lowers the resolution as you zoom out? That would be perfect for me, my data is large but you only need the high resolution once you're zoomed in to a specific area.


texts: Fixed distance between two contours - Qgis 2.14.1

texts:
There is a possibility, when you select a contour line and give a fixed length of line, this line connects automatically to the next contour lines? I have two other questions on this subject but I do not have an adequate response.
This is very important for the design of forest roads, to select an optimum zero line, depending on the longitudinal inclination of the forest road. (1-12%). Qgis 2.14.1 Essen

Connect two contour lines with calculated line length

Connect two points with fixed longitudinal slope over DEM - QGIS


texts:
How can I find a max cell value from a stacked raster.

Rmax <- maxValue(RAD1998.all[[1]])

works fine but

Rmax <- maxValue(RAD1998.all[[2]])

gives NA.

Of course not in a stacked raster.

Here is my code:

RAD1998 <- raster(paste(getwd(), "/1998bil/1998ASC5min_ppt_spas1214_0001_19980202_0810_UTC.asc.bil", sep = ""))
list.ras <- mixedsort(list.files(paste(getwd(), "/1998bil/", sep = ""), full.names = T, pattern = ".asc.bil"))
RAD1998.all <- stack(list.ras)


texts: Renaming layer with chart pie symbology in TOC of ArcMap?

texts:
I work with ArcGIS 10.3 and when i try rename layers point1 point2 in the TOC but i can't. I using chart pie symbology:

This is the symbology tab:


texts: Spatial joining two polygon objects in R?

texts:
I want to combine voter precinct data with school zone data. Precinct boundaries do not line up neatly with school zone boundaries in all cases.
Code look something like this:
library(rgdal)
library(tigris)
library(sf)

precincts <- readOGR(
  dsn = "input/SOEShapefiles2020/SOEShapefiles2020.shp",
  layer = "SOEShapefiles2020"
)

precinctDemos <- read.csv("input/data/Party registration in voting precincts.csv")
mPrecDemo <- geo_join(
  spatial_data = precincts,
  data_frame = precinctDemos,
  by_sp = "PRECINCT",
  by_df = "Precinct",
  how = "inner"
)

elemZones <- readOGR(
  dsn = "input/2020-21_Attendance_Boundaries/2020-21_ElemAttendance.shp",
  layer = "2020-21_ElemAttendance"
)
elemZones$MSID_ELEM <- as.numeric(elemZones$MSID_ELEM)
attendance <- read.csv("elementary school attendance.csv")

mElem <- geo_join(
  spatial_data = elemZones,
  data_frame = attendance,
  by_sp = "MSID_ELEM",
  by_df = "School.",
  how = "inner"
)

So mElem and mPrecDemo are each "Large SpatialPolygonsDataFrame" objects. How do I perform a spatial join to bring the two objects together? I want to be able to estimate the number of voters in each school zone.


texts: Pan & Change Zoom to show Fusion Tables Query Results in Google Maps JavaScript API V3

texts:
I've created this map --  https://s3.amazonaws.com/pdb.files/lookup_auto2.html -- using Google Fusion Tables and Google Maps API. I'm attempting to modify the code so that the map will pan and  zoom after executing the query. It seems that I need to use either the fitBounds function to display the bounds of the results, or something simpler like this:

// Change the center and zoom of the map
      map.setCenter(results[0].geometry.location);
      map.setZoom(10);

I'd prefer the second option so that I can set a consistent zoom level for each query.  Any pointers or references for completing this are much appreciated.


texts: Using shared database authentication settings for organization environment?

texts:
My company is willing to share common QGIS3 configuration parameters between its users (ie. proxy configuration, WFS/WMS connections definitions, database connection definitions).

For doing this, the documentation explains how to use a global settings file (ie. a customized qgis_global_settings.ini file), which works perfectly fine. For the case of database connections (Postgis & Oracle in our case), I configured the connections this way:

[PostgreSQL]
connections\myconnection\service=
connections\myconnection\host=my_company_database_host
connections\myconnection\port=5432
connections\myconnection\database=my_database
connections\myconnection\username=xxxxxxx
connections\myconnection\password=yyyyyyyy
connections\myconnection\authcfg=
connections\myconnection\saveUsername=true
connections\myconnection\savePassword=true
(...)

This raises 2 problems:

'username' and 'password' are plaintext/unencrypted.
for this to work, I had to switch 'saveUsername' and 'savePassword' to TRUE, because we want to avoid our users the burden of entering the credentials by themselves.

I believe that saveUsername and savePassword will save the credentials in any project file using this database connection. This is a problem if the project files are sent outside the company.

Another approach:

QGIS3 provides an authentication manager. It stores encrypted passwords fields, in a sqlite database file (qgis-auth.db). My understanding is that, in some way, this auth database can be linked with the above configuration file using the authcfg field.

Is there any way to use the qgis-auth.db in an organization environnement ?

For example,

sharing a qgis-auth.db file on network
using the Python API to fill a local qgis-auth.db with default credentials when deploying QGIS for the first time


texts: PyQGIS error: proj_create_from_database: Cannot find proj.db

texts:
I have problem with QGIS. When I run my code in QGIS Python console it works but when I try in console I'm getting an error. Code:
import sys

from qgis.core import QgsApplication, QgsProcessingFeedback
from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis', True)
qgs = QgsApplication([], True)
qgs.initQgis()

# Add the path to processing so we can import it next
sys.path.append(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis/python/plugins')

import processing
from processing.core.Processing import Processing

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
feedback = QgsProcessingFeedback()

processing.run("native:joinattributesbylocation", {'INPUT':'C:/Users/ja/Inne/Desktop/bdot10k_bubd_a/bdot10k_bubd_a.gpkg','JOIN':'C:\\Users\\ja\\Inne\\Desktop\\POLSKA\\PRG_PunktyAdresowe_POLSKA.shp','PREDICATE':[1,5],'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],'METHOD':1,'DISCARD_NONMATCHING':True,'PREFIX':'','OUTPUT':'C:/Users/ja/Inne/Desktop/projekt2.gpkg'})

What is the problem?
Error: >>>
... from qgis.core import QgsApplication, QgsProcessingFeedback
>>> from qgis.analysis import QgsNativeAlgorithms
>>>
... QgsApplication.setPrefixPath(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis', True)
>>> qgs = QgsApplication([], True)
>>> qgs.initQgis()
>>>
... # Add the path to processing so we can import it next
... sys.path.append(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis/python/plugins')
>>>
>>> import processing
>>> from processing.core.Processing import Processing
>>>
... Processing.initialize()
>>> QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
proj_create_from_database: Cannot find proj.db
True
>>> feedback = QgsProcessingFeedback()
>>>
>>> processing.run("native:joinattributesbylocation", {'INPUT':'C:/Users/ja/Inne/Desktop/bdot10k_bubd_a/bdot10k_bubd_a.gpkg','JOIN':'C:\\Users\\ja\\Inne\\Desktop\\POLSKA\\PRG_PunktyAdresowe_POLSKA.shp','PREDICATE':[1,5],'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],'METHOD':1,'DISCARD_NONMATCHING':True,'PREFIX':'','OUTPUT':'C:/Users/ja/Inne/Desktop/projekt2.gpkg'})
proj_create_from_wkt: Cannot find proj.db
proj_create_from_database: Cannot find proj.db
ERROR 1: PROJ: proj_identify: Cannot find proj.db
proj_create_from_wkt: Cannot find proj.db
proj_identify: Cannot find proj.db
proj_get_authorities_from_database: Cannot find proj.db
ERROR 1: PROJ: proj_identify: Cannot find proj.db


texts: Get Parameter as Text and select layer by attribute

texts:
I have written a python script tool wherein a user specifies a feature class (e.g., Landuse) containing different attributes (e.g., parcel size, landcover etc.), name of the attribute field from which features will be selected (e.g., "landcover") and the name of the attribute to be selected (e.g., "forest" or "buildings"). The end result should be that the script selects all the features with the specified attribute and create a new feature class out of it. For this I have 3 input parameters, as shown below, and the script receives the input with GetParameterAsText().

My code block looks like follows:

FClass = arcpy.GetParameterAsText(0)
Field = arcpy.GetParameterAsText(1)
Feature = arcpy.GetParameterAsText(2)

arcpy.MakeFeatureLayer_management(FClass, "FclassLayer")
arcpy.SelectLayerByAttribute_management("FclassLayer","NEW_SELECTION","'Field' = 'Feature'")

arcpy.CopyFeatures_management("FclassLayer", "Feature1.shp")

When I run the script tool, it doesn't produce any error message and creates Feature1.shp, but the shapefile is simply a copy of the specified feature class. That means, SelectLayerByAttribute is simply selecting the whole feature class irrespective of the specified attribute.

Being a beginner with python scripting, I am not being able to locate what am I doing wrong here. Furthermore, I tested the same script with hard coded values instead of getting it from user with GetParameterAsText, and there it runs perfectly and gets me the output as expected.


texts: Update value in domain and same should reflect in attribute field

texts:
I have a SDE GDB, where i have no. of domains created. we have created domain for state names, which is used by all feature classes to update their respective state name.

Now i have discovered that spelling / name of some of the states are not correct. if i correct these spellings or names in domain. will it reflect in the attribute of all feature classes or is there any other way forward?

There are many thousand records in the GDB.


texts: Cannot save GPX files in ArcGIS Online and make them to perdure in My Content

texts:
I am trying to use a sample.gpx containing 2,913 records and lat, lon, name, cmt and type per record, but uploading it into ArcGIS Online (Public account) and works fine at the very beginning, in essence I can see my points displayed in the map and/or change the symbology and so on, but after Saving or Signing Out -still I haven't isolated the cause- when I move back to the map for editing, I cannot see the points displayed any more! no matter how long I wait for the map to load (in case it is a server response issue) and I have to re-upload them again as a new layer, otherwise I don't know how to bring them back.

It is very frustrating starting over and over again, any ideas?


texts:
I downloaded a HDF file from here (MOD13Q1). Upon viewing the content I found that the file has 4800X4800 rows/columns(Resolution 250 meters). Now I tried to view the same file in QGIS. Upon zooming in I saw rectangular polygons having the EVI values. But in the file, I didn't saw any geometry column.

Can anyone explain how QGIS interpereted the rectangular polygon(pixel) ?


texts: Unable to display FeatureLayer using ArcGIS server

texts:
I can't figure out why this FeatureLayer doesn't show up on top of my base map. I'm using ArcGIS Server 10.0 and the ArcGIS JS API 2.6.

Here's the information from the service's REST parameters page:

Feature Layer: Inspection Locations (ID: 5)
Display Field: NAME
Type: Feature Layer
Geometry Type: esriGeometryPoint

Min. Scale: 250000
Max. Scale: 0

Extent:
XMin: -8879044.1973
YMin: 4938432.3627
XMax: -7999181.0204
YMax: 5622923.0498
Spatial Reference: 102100 (3857)
.
.
.
Supported Interfaces:   REST
Supported Operations:   Query   Add Features   Update Features   Delete Features   Apply

Edits

...and here's the code:

  <script type="text/javascript"> var djConfig = { parseOnLoad: true }; </script>
  <script type="text/javascript" src="http://serverapi.arcgisonline.com/jsapi/arcgis/?v=2.6"></script>
  <script type="text/javascript">djConfig = { parseOnLoad:true };</script>

  <script type="text/javascript">

  dojo.require("esri.map");
  dojo.require("dijit.layout.ContentPane");
  dojo.require("esri.layers.FeatureLayer");
  var map;

  function init() {

    var initExtent = new esri.geometry.Extent({
      "xmin": -8879044.1973,
      "ymin": 4938432.3627,
      "xmax": -7999181.0204,
      "ymax": 5622923.0498,
      "spatialReference": {
        "wkid": 102100
      }

    });
    map = new esri.Map("map", {extent: initExtent });

    var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer");

    var featureLayer = new esri.layers.FeatureLayer("http://gis4/ArcGIS/rest/services/NY/NYeFind/FeatureServer/5",{
     mode: esri.layers.FeatureLayer.MODE_SNAPSHOT,
     outFields: ["*"]
   });

    map.addLayers([basemap, featureLayer]);

    var resizeTimer;
    dojo.connect(map, 'onLoad', function(theMap) {
      dojo.connect(dijit.byId('map'), 'resize', function() {
        clearTimeout(resizeTimer);
        resizeTimer = setTimeout(function() {
          map.resize();
          map.reposition();
        }, 500);
      });
    });

    } // end init

    dojo.addOnLoad(init);

Oddly enough, if I hit the "View in ArcGIS.com Map" link, it works perfectly.


texts:
I tried to plot Variogram by this code from gstat packages:

v <- variogram(log(zinc) ~ 1, meuse)
plot(v)

and here is the output:

But I want to add line to point (draw variogram line).

How I can add line to the plot?


texts:
I'm working on a project where I need to snap an existing line layer to a different line layer. These lines represent streams and rivers. I'm working within QGIS Maderia 3.4 LTR.

I've tried using the "Snap Geometries to layer" tool but it doesn't appear to snap the input layer to the reference layer, the input layer continues to stay in the same position, no matter the tolerance I put in (10ft, 15ft, 100ft, 500ft).

I've figured out that I can use the vertex tool, enable snapping and manually do a lot of it. With the snapping and the vertex tool, I've played around with the type also (vertex, vertex and segment, segment).

With the manual snapping, I've been having particular trouble with the confluences, where one stream merges with another. This is an area of the snapping that especially needs to be done precisely. If this process does need to happen manually, any hints for that particular task?

The dataset is large and I was hoping to find a less tedious process. Especially since I need to be able to snap or create a snapped layer to within a 15ft precision level (currently with the manual snaps, I've been using a buffer as a guide).


texts: How to calculate break point values and ranges for point density heatmap?

texts: I have some questions about the point density analysis in ArcGIS. I am interested to know the concept and formula about how to come out with all the break point values and ranges of the different colour bands for point density heatmap. I hope any of you can provide me a rough idea about the concept and the formula to use. I want to know the rough Value/Range I can get using the formula.


texts: Euclidean distance units for WGS84 files?

texts: I have a lake raster file and I want to calculate the euclidean distance. My raster file is in EPSG 4326 and I use the r.grow.distance from GRASS GIS through QGIS plugins. What are the units from the r.grow.distance for a file with EPSG 4326 and how can I transform the distance to meters? The area is between 50N to 60N.


texts: Rasterio: select subset of georeferenced image by lop left coordinates and height, width

texts:
I want to subset a 10m sentinel dataset using rasterio window modul.
I have tried using the window = from_bounds(...) part in .read() however it requires for me to put in the left top right bottom coordinates.
I want to subset using only the top left coordinates and the height and width of subset in pixels.
Here is what I have tried so far
import rasterio
from rasterio import plot # for plotting datasets, may not be of use
from rasterio.plot import show_hist # for showing histogram distribution of bands in a daatset, may not be of use
from rasterio.windows import Window , from_bounds, crop# for clipping out only a window from dataset use in part 2
import matplotlib.pyplot as plt #for plotting images
import numpy as np # for doing computations
import os  # to read file paths to the image data
from sklearn import cluster
%matplotlib inline

path_10 =  R'./GRANULE/L2A_T31UET_A025788_20200530T105134/IMG_DATA/R10m'
l_10 = os.listdir(path_10)
tostack_10 = l_10[1:5] #take bands B G R NIR for 10m data
stackpath_10 = [os.path.join(path_10,i) for i in tostack_10]
ds_10 = [rasterio.open( stackpath_10[i], driver='JP2OpenJPEG' ) for i in  range(len( stackpath_10 ))] # list of datasets

bands_10= [] #array holding all bands as numpy arrays

# window_= Window( 601200.0, 5773695.0, 700, 500)

for i in range( len(ds_10)):
    bands_10.append( np.array(ds_10[i].read(1 ) , dtype='int32' ) )
    ds_10[i].close()

bands_10 = np.array(bands_10 , dtype='int32') # typecasting from unsigned int to signed
bands_10.shape

This gives me a shape of (4, 10980,10980), which is teh shape of the total scene and not the subset


texts: Organize all dialogs in one package QGIS 3 Plugin

texts:
To organize better the code in my QGIS 3 Plugin, I tried to put all dialogs in one package called dialogs, but when I open QGIS it says ModuleNotFoundError: No module named 'dialogs'

This is the code I used to import the Dialog class in the "main" file, where there's run method (luniPlugin.py in this case)

from dialogs.select_feature_dialog import SelectFeatureDialog

(and the same for the others dialogs)

This is the folder where the project is:

And, finally, this is the content of the folder "dialogs":

Since python compiler has not given me errors, I'm wondering if it's possible to do this "unpacking" in a QGIS Plugin or if it's simple my mistake.


texts: Align Raster with boundary shapefile: selecting right proj4 R

texts:
library(maps)
library(mapproj)
library(mapdata)
library(rgeos)
library(maptools)
library(sp)
library(raster)
library(rgdal)

Drawing inspiration from here https://stackoverflow.com/questions/10763421/r-creating-a-map-of-selected-canadian-provinces-and-u-s-states, I would like to align the raster and shapefile below. Information from the data providers say datis on a polar stereographicprojection. However, I could not get it to my shapefile to align with the rasterLayer. Here is my attempt so far:

dat file~1 KB

#saveRDS(dat1,'dat.rds')

dat<- readRDS('dat.rds')

# can0<-getData('GADM', country="CAN", level=0) # Canada
can1<-getData('GADM', country="CAN", level=1) # provinces

## Specify a geographic extent for the map
## by defining the top-left and bottom-right geographic coordinates
mapExtent <- rbind(c(-156, 80), c(-68, 19))

## Specify the required projection using a proj4 string
## Use http://www.spatialreference.org/ to find the required string
## Polyconic for North America
newProj <- CRS("+proj=poly +lat_0=0 +lon_0=-100 +x_0=0
            +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

## Project the map extent (first need to specify that it is longlat)
mapExtentPr <- spTransform(SpatialPoints(mapExtent,
                                         proj4string=CRS("+proj=longlat")),
                           newProj)

## Project other layers
can1Pr <- spTransform(can1, newProj)
#us1Pr <- spTransform(us1, newProj)

## Plot each projected layer, beginning with the projected extent
extent(dat)=extent(can1Pr)
plot(mapExtentPr, pch=NA)
plot(dat, add=TRUE)
plot(can1Pr, border="black", col=NA, add=TRUE)

#===================

Additional information about the data:

The CANGRD grid is in polar stereographic projection with a 50 km spatial resolution. The grid is a 125 (columns) by 95 (rows) matrix, where the SW corner (0,0) is at 40.0451°N latitude and 129.8530°W longitude. The projection is true at 60.0°N and centered on 110.0°W. The file ‘CANGRD_points_LL.txt’ lists the latitudes and longitudes for each grid point.

The general format of the ‘YYYYDD.grd’ file is:
Id – ‘DSAA’ identifies the file as an ASCII grid file
nx ny - nx is the integer number of grid lines along the X axis (columns)
        ny is the integer number of grid lines along the Y axis (rows)
xlo xhi - xlo is the minimum X value of the grid
          xhi is the maximum X value of the grid
ylo yhi - ylo is the minimum Y value of the grid
          yhi is the maximum Y value of the grid
zlo zhi - are the minimum and maximum Z values of the grid.

data source:
ftp://ccrp.tor.ec.gc.ca/pub/EC_data/CANGRD/


texts: Snapping points to lines in QGIS/PostGIS?

texts:
I'm attempting to snap approximately 900 points to a line table.

I've tried recreating what I have found from a couple of PostGIS answers but not having much luck. To test, I've created two simplified tables; one for points (ptest) and one for the lines (ltest). The test tables only consist of an ID column.

CREATE TABLE line_attribute_point AS
    SELECT
        distinct on (ltest.id) ltest.*,
        ST_ClosestPoint(ltest.geom, ptest.geom) as snapped_point
    FROM
        ltest,ptest
    ORDER BY ST_Distance(ltest.geom, ptest.geom);

The main error seems to be the ORDER BY.

Any ideas?


texts: Is accessing an mxd from an Msd based SOE safe?

texts:
I have a server object extension (SOE) that returns a list of legend items based on symbols that fall within a specified extent.  This allows the client to update a legend to reflect on the symbols that fall within the map extent.

To do this requires access to renderers.  This is simple with with mxd based mapservices, however Msd based mapservices don't support access to renderers.  As a workaround I'm  using IMapDocument to open the mxd and keep a reference to the layout as a member variable.

Is this safe?

Support for mxd mapservices will be deprecated at 10.1, so I really need to figure this out.

Update

This video from the Dev Summit doesn't explicitly address this question. However, given how the printing gp service uses mxd files, this implies mxd access is safe.  If it isn't, I suppose I could re-deploy my SOE as a gp service.

Update2

Here's a related question on CodeReview.Stackexchange.


texts: How to integrate and plot both heatmap and quivers/arrows for different timestamps using folium?

texts: I want to plot wind velocity heatmap and wind direction quivers together for different timestamps on top of a map using folium. I could not find any plugin for this integrated operation. Can you please help me sort this out?


texts: How would I go about re-projecting the Blue Marble raster maps into a robinson projection using free tools?

texts: I would like to convert the NASA Blue Marble maps into a robinson projection.


texts: Creating DTM for two surfaces with point data and calculating volume between them using ArcGIS Desktop?

texts:
I have point data for a ground surface level which is captured in a grid of 35m spacings.  I also have point data for the new surface level after soil has been removed which is captured at random locations using a mobile GPS.  I want to create a DTM for both surfaces and then calculate a volume for the amount of soil removed.

Is this possible with ArcGIS Desktop 10.0?


texts: ArcGIS 10.2 layers in data view but not in layout view or table of contents

texts:
I have spent considerable time editing line feature classes in one of my file geodatabases in an ArcMap 10.2 map file. Today, I returned to it, opening the map file, and none of the layers appeared in the contents panel (list by drawing order). I checked out the contents tabs -- nothing shows in 'list by source', but I can finally see my database objects in 'list by visibility' and 'list by selection.' I can't open their properties, attribute tables, etc. I could still see the edits that I made there, too, so that's good at least.

Switching to layout view, nothing was there. Back to data view. Still there. Layout view, nothing. I tried toggling drafting mode. No good. Coordinates? Should be there.

I decided to move right ahead with topology to make sure my lines were all lined up properly so I could convert my lines to polygons. I enter edit mode, say I want to edit my geologic lines (weird that it shows up there, but not in the contents), open up my topology feature, select the features I want to check, and I can't open the error inspector.

Fine, so something's off and it looks like I can't just ignore it. I'm used to Arc throwing errors.

I opened ArcCatalog and opened my map document through that. Everything shows up fine in contents, etc but all of the work I did the other day is GONE!

I close all of Arc and open up the map document directly like I did earlier. My edits are there, but I have the same problems as before. I load my old line featureclass into the map document. It only shows up in the contents. I switch to layout view, and I can see the lines. But wait... where are my edits!? Switch back to data view. Ok, there they are. Back to layout view. Edits not there anymore.

I realize that I'm rambling, but I'm not quite sure how to communicate this issue. How do I recover my edits? They've clearly been saved somewhere, but not in the file geodatabase where they should have been saved... or are they? I can't check the layer properties to make sure that the data source is correct, and I'm more than half-convinced that what I see is only a "ghost" of my data, but really don't want to spend another day drawing all of those lines again.


texts: Getting geometries of hurricane's 'cone of uncertainty' using shapely?

texts:
I'm trying to plot hurricane forecast maps using python. I have several forecast positions derived from official advisories, and interpolate them into a smoothed curve, then draw a polygon of 'cone of uncertainty' based on the curve. Example:

Basically, 'cone of uncertainty' is the footprint of a moving and enlarging circle. I've tried many approaches but none of them is good enough. My current approach is to produce ~100 circles based on the interpolated curve, and make a compound polygon using cascaded_union method in shapely.

import numpy as np
from shapely.geometry import MultiPolygon
from shapely.ops import cascaded_union
from scipy.interpolate import interp1d
# x, y: coords of forecast position
y = [18.3, 19.2, 20.0, 20.4, 20.7, 21.3, 21.6, 21.5, 20.8, 20.8, 21.5]
x = [111.3, 111.2, 110.9, 110.5, 110.2, 110.5, 110.0, 109.2, 109.4, 110.3, 111.8]
# r: radius of uncertainty
r = [0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.5, 0.5, 0.5]
hours = [0, 6, 12, 18, 24, 36, 48, 60, 72, 96, 120]
# interpolate
points_num = 100
interp_hours = np.linspace(min(hours), max(hours), points_num)
x = interp1d(hours, x, kind='cubic')(interp_hours)
y = interp1d(hours, y, kind='cubic')(interp_hours)
r = interp1d(hours, r, kind='linear')(interp_hours)
# make polygon
thetas = np.linspace(0, 2 * np.pi, 360)
polygon_x = x[:,None] + r[:,None] * np.sin(thetas)
polygon_y = y[:,None] + r[:,None] * np.cos(thetas)
polygons = MultiPolygon([Polygon(i) for i in np.dstack((polygon_x, polygon_y))])
polygons = cascaded_union(polygons).buffer(0)

But it looks nasty near the starting point:

Increasing the number of circles can only partly address the problem and it takes more time. So I wonder if there is a beautiful, efficient and pythonic way to make the 'cone of uncertainty'? Note that hurricanes may change its directions abruptly, and even be stationary!


texts: QGIS Grass Plugin fails with r.cost.vect

texts:
I tried the following in QGIS with the Grass Plugin:

r.cost input=preCostSurface@albi start_points=OSMroads@albi output=cost0109 max_cost=0

And I get the following error:

Dev note: Adapted sites library used for vector points. (module should be updated to GRASS 6 vector library)

￼No start points
Finished with error

What is causing that and how can I fix it?

Edit: It seams like this is a general Grass 6.4.3 issue and will be resolved in Grass 7. Does anybody know a workaround for that problem meanwhile?


texts: Import CSV column into existing feature using arcpy

texts:
My goal is to take a CSV that contains 50 rows and one column and import the data into a feature class with over 3,000 rows using arcpy.  I want to pick a random row from the CSV and then iterate through the rows in the featureclass.

I will need to use this workflow constantly, but am unsure of how to read the CSV into arcpy.


texts: GPS with QGIS--what happened to the capture point?

texts:
I have been using QGIS for years but when v. 2 came out and important functionality seems to no longer work.  In 1.8 the "capture point" or "capture line" buttons on the GPS panel would work.  Ever since 2.0 I have been unable to make these work.  What happened?  Maybe it is something different than QGIS because I also have switched to using windows 8 tablets in the field environment where I use the system.

btw--a related question, QGIS doesn't seem to be able to capture Z information in any case.  Capturing Z for points, to me, is an essential feature that seems missing.  Maybe I'm missing something here that is in v. 2 that I just don't understand.
TP


texts: How to compute buffers for certain features?

texts:
I used the following snippet of codes for finding buffer.

def buffer(self):
    layer = self.iface.activeLayer()
    distance = float( self.dlg.DistancelineEdit.text() )
    city = self.dlg.Citycombo.currentText()
    expr = QgsExpression('"name_1" = \'' + city + '\'' )
    it = layer.getFeatures( QgsFeatureRequest( expr ) )
    feature = it.next()

    # Perform a buffer on the feature geometry
    myBufferPolygon = feature.geometry().buffer(dist, 50)

It is raising two errors

feature = it.next()
    StopIteration.

2.distance = float( self.dlg.DistancelineEdit.text() )
ValueError: could not convert string to float:


texts: Why does Add Field to Layer File give "ERROR 000840: The value is not a Feature Layer"?

texts:
I'm trying to add a field to a layer file, but it always gives an error. What I did was making a closest facility layer and save the result as a layer and after that I tried to add a field, but then it gives an error. It says it cannot add a field if it isn't a raster layer, ...

I tried to avoid this by making a shapefile (manually) and this can work, but is there a method to do this in Python? So convert a route from the closest facility layer to a shapefile using Python? It always says:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Feature Layer.


texts: Overlaying ArcGIS REST API on Google Maps?

texts:
For example I have ArcGIS REST API as
us high way state city county

https://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Specialty/ESRI_StateCityHighway_USA/MapServer

I want to see it on Google Maps.

How do I do it?


texts: Reducing points but maintaining attributes from original shapefile

texts:
I have a shapefile with points but there are a lot of points and I want to reduce the amount.
I have tried "Points thinning" which automatically reduces the amount of points but it loses all the attributes/features/values from the original file (which I want to keep). In the 'Thinned Points' there are columns: 'Count', 'Mean', 'Minimum', 'Maximum', 'StdDE' which I don't want.
How can I reduce points but maintain the attributes from the original file?


texts:
What would be the steps or process to migrate data from an ArcFM system to a vanilla installation of ESRI ArcGis?
We tried loading the SDE data into an ArcGIS installation without ArcFM and received an unable to create object class extension error.


texts: GeoDjango admin: how to use a WMS background other than the default OSM or VMAP0

texts:
The standard GeoDjango admin has a basic VMAP0 background. OsmGeoAdmin adds an OSM background.

How do I replace the background map with another WMS?

In the source code of the admin OpenLayers template there are {{ wms_url }} and {{ wms_name }} variables but I cannot find where these can be specified to customize it.

Or should I just completely override the default admin template?


texts: Return field names using Arcade function in ArcGIS Online

texts:
I'm trying to write an expression for a pop up in ArcGIS Online that will display the name of the fields if there is a specific value under those fields for each feature. I adapted the following script that I found online.

var i = 0
var features = [];

function addvalue(feat) {
    if (!IsEmpty(feat)) {
        features[i++] = feat;
    }
}

addvalue($feature.schoolCommGardenQuarter);
addvalue($feature.smarterLunchQUARTER);
addvalue($feature.schoolWellnessQuarter);
addvalue($feature.keysNAPSACCquarter);
addvalue($feature.healthConcessionsQuarter);
addvalue($feature.foodPantryBackpackQuarter);
addvalue($feature.farmMarketorPopUpQuarter);
addvalue($feature.rydQUARTER);
addvalue($feature.directEd);
addvalue($feature.onetimeRYDEventQuarter);

return Concatenate(features, ";")

This runs without issue but it returns the values as a concatenated list (such as [YES;NO;YES;YES]). I want to edit the above script so that it returns a text string of the field names that have YES as a value.

I have since edited my script so that the if statement evaluates for the 'yes' value. I still am having trouble figuring out how to access and return the field names rather than the 'yes' values.

var i = 0
var features = [];

function addvalue(feat) {
    if (feat == 'Yes') {
        features [i++] = feat;
    }
}

addvalue($feature.schoolCommGardenQuarter);
addvalue($feature.smarterLunchQUARTER);
addvalue($feature.schoolWellnessQuarter);
addvalue($feature.keysNAPSACCquarter);
addvalue($feature.healthConcessionsQuarter);
addvalue($feature.foodPantryBackpackQuarter);
addvalue($feature.farmMarketorPopUpQuarter);
addvalue($feature.rydQUARTER);
addvalue($feature.directEd);
addvalue($feature.onetimeRYDEventQuarter);

return Concatenate(features, ";")


texts: Cesium js, Draw line binding a label to a position

texts:
I would like to draw a line binding a entity to its label with an offset.
CesiumJS allows to offset the label, however its not possible to draw a line (or polyline) from a position to an offset like the red line in this image.

How can i do it? any sugestion?

i'm using pixel offset. but there is no problem to use eye offset

 labels.add({
        position: Cesium.Cartesian3.fromDegrees(-75.1641667, 29.9522222),
        text: 'Another label',
        pixelOffset: new Cesium.Cartesian2(100,-100)

    });


texts: Formating PDAL Pipeline JSON for crop filter


texts: How can I change layer names in my legend when using layers with representations? Arcgis 10.3.1

texts: When using a Feature Class layer with symbology converted to a representation there is no option in the symbology tab of the layer's properties to enter legend text. I can't see any options in legend properties. Does anyone know how to do this?


texts: Select by Attributes - Strings with < 5 characters

texts:
I have a column of zip codes which are strings.  Most of them are the usual 5 digit zip codes, but some are only 3 digits or missing a digit. (keep in mind, I need this column as strings for other reasons)

ex.

  ZIP
  234
  4566
  12345

I basically want to select everything that has less than 5 characters for further investigation.  I'm sure there's an easy way to do this, but I'm having a difficult time figuring it out.


texts: Seeking open source software package for Remote Sensing?

texts:
Here are many questions with great answers about open source GIS software.

I am wondering, what is the best open source software package for Remote Sensing? I would like to learn it and to use in my work.

I used to work with IDRISI, and I've heard about Erdas and ENVI, but they all are not free. I am looking for a free and powerful leader, like Qgis for GIS or R for statistics. With powerful classification, segmentation, Fourier, filters, PCA, etc.

Can anyone please advise me a good free RS software? What are the capabilities, user friendly or with command line? Do any comparison matrices exist?


texts: How can I assign Raster2 pixel value to Raster1 where they overlay?

texts:
I want to assign Raster2 pixel value to Raster1 (bigger) where they overlay. Both of them have the same cell size.

How can I accomplish this process in ArcGIS Desktop? (Like Raster Calculator)


texts: Taking large netcdf and converting it to interpolated values in different set of coordinates?

texts:
I've got a large netcdf (a small portion of the data is here) with the data in an even grid associated with an Albers Equal Area Conic projection.  I've got the corresponding latitudes and longitudes here and here.

I want to take the information in this netcdf, and get interpolated values along a regular lat/lon grid, from 18.2W to 53E and 32.8S to 20N.  For output, I'd like a 3-dimensional array with the x, y, and z axes being lon, lat, and time.  Or a set of rasters that I can coerce to a 3d array in R.

How can this be done?

I know that I could write an R script where for each point in my regular lat/lon grid I identify the four nearest surrounding Albers points and take their inverse-euclidian-distance-weighted average.  But this will be painfully slow to compute and difficult to program.

Are there any programs or tools that can help me to batch-interpolate a netcdf with a deep time dimension?

I'll have >100 of these to do, and each grid is 1152x1152.

I seek solutions in R, or anything that can be run in some sort of batch mode.

This was a case where I simply didn't know the names of the tools to use.  Apparently, the "projectRaster" function in the raster package (which relies on gdal) does what I need once I coerce my netcdf slices into raster objects.  I had to do a bit of digging in metadata and obscure documentation to figure out proj4 syntax, but that was the extent of the difficulties.


texts: Why are steps out of order when running model in Batch?

texts:
I have an issue using an ArcMap model in batch mode. The purpose of my model is to obtain the result of an OD cost matrix from a selected network and a set of origins and destinations. My goal is to be able to get these results without having to add the data to the display in ArcMap (time and resource consuming), but most of all, I want to be able to automatically process several OD cost matrix one after another. For these reasons, I rebuilt the OD cost matrix routine using ArcMap's Model Builder.

If I test the model on only only network, it works great. However, when I try to batch the process where each time I have a different network with different locations, I have a problem since ArcCatalog seems to mix the different steps. For example, if I batch with only two different networks...

Instead of having:

1. Create OD cost layer (1)
2. Load Origins 1 in network 1
3. Load Destinations 1 in network 1
4. Solve OD cost matrix (1)
5. Create OD cost layer (2)
6. Load Origins 2 in network 2
7. Load Destinations 2 in network 2
8. Solve OD cost matrix (2)

ArcCatalog does this:

1. Create OD cost layer (1)
2. Create OD cost layer (2)
3. Load origins 1 in network 2
4. Load origins 2 in network 2
5. Load destinations 1 in network 2
6. Load destinations 2 in network 2
7. Solve OD cost matrix (1)
8. Solve OD cost matrix (2)

I thought that when batching a model ArcCatalog would do the first task completely before even looking at the second. Why are the steps of my model mixed up?

Thanks for your help!


texts:
I am trying to do kriging by gstat package through R.

Here is my Code:

krg<-krige(var ~ 1,data,grid,model=vgm,nmax=30)

But each time I recive this warning:

1: In predict.gstat(g, newdata = newdata, block = block,  ... :

 Warning: Covariance matrix (nearly) singular at location [1,4501,0]: skipping...

and when I try to plot it by spplot this error appears:

Error in seq.default(zrng[1], zrng[2], length.out = cuts + 2) :

'from' cannot be NA, NaN or infinite
In addition: Warning messages:
1: In asMethod(object) :
  complete map seems to be NA's -- no selection was made
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf

Here is the full code:
You can download my data through this link: Mydata
and here is my full code:

library(gstat)
library(sp)
library(raster)
coordinates(Mydata) <- ~East + North
proj4string(Mydata) <- CRS("+proj=utm +zone=10+datum=WGS84")
ext <- as(extent(Mydata), "SpatialPolygons")
r <- rasterToPoints(raster(ext), spatial = TRUE)
proj4string(r) <- proj4string(Mydata)
v1 <- variogram(cpue_female ~ 1, Mydata)
vgm_mod_omni <- vgm(psill = 0.0003, model = "Mat", range = 50000, nugget = 33, kappa = 0.5)
G1 <- krige(cpue_female ~ 1, Mydata, r, vgm_mod_omni, nmax = 30)
gridded(G1) <- TRUE
spplot(G1)


texts:
I want to complete a code that should create same OAuth2 authentication on many other computers but with different usernames. I already tested creating basic authentication and it works fine.

am  = QgsApplication.authManager()
cfg = QgsAuthMethodConfig()
cfg.setId('test')
cfg.setName('test')
cfg.setMethod('Basic')
cfg.setConfig('username', "username")
cfg.setConfig('password', "password")
am.storeAuthenticationConfig(cfg)

But when I wanted to add OAuth2 authentication I stuck at adding bunch of parameters. Need to point Resource Owner grant flow and specify url, cliend id, client secret, username and password.

Like here:

Also tried to read configs of existing authentications to puth them into new configs but everytime got an empty configMap.

am  = QgsApplication.authManager()
for name, conf in am.availableAuthMethodConfigs().items():
    print(name)
    print(conf.configMap())

Any thoughts?


texts: Using Iterator with ModelBuilder for clipping and splitting in ArcMap?

texts:
My input feature class is a county-wide dataset and the input clip feature is a township’s ‘section.’  I have 576 different ‘sections’ in which to clip the county-wide data out to.  Though I dont have 576 different section feature classes to clip by, I was just using a definition query to select which particular section to clip out to.  I would like to use ModelBuilder to streamline clipping 576 sections of countywide dataset.

How would I set up ModelBuilder to incorporate an iterator into the model for this particular case?

The above model successfully split a township's contours into 36 sections. I need to run it 15 more times to achieve my end goal. I need to add an iterator to batch split out the remaining 15 townships. I put an orange background on the items that need to have the expression or output filename changed when a different township is being run through the model.


texts: Understanding raster custom transparency options (making null data 100% transparent) in QGIS?

texts:
I have a raster with a large null area (see the solid black area in the image below). The "Band 1" value throughout this null area is -3.402823e+38.

I'm trying to make sense of the result of my setting the "custom transparency options" so that "Band 1" values from 0 to 0 are 100% transparent (see settings below). Please note that at no point am I making changes to the default Style properties of the raster.

As a result of these settings, the null area does become transparent but the rest of the raster now has strange banding (see result below). My expectation was that the only visual difference between the initial raster and the result raster would be 100% transparency in the null area.

Am I missing something?

EDIT: Per request, here are the the Style settings:

And I have set '0' as an Additional no data value, as shown below.

I have further discovered that I achieve the same strange banded raster result even if I delete the custom 0 to 0, 100% transparency option. To achieve this strange banded raster, I only need to select "Band 1" as the transparency band as shown below.


texts: Python script causing python.exe crash and wrapped C/C++ object failure

texts:
I am working on automating some of the normal processes that I want to do and I have run into some problems I have not found solutions for after a couple days of searching.

I am relatively new to python / QGIS.

QGIS version 2.16.3 (still occurs in 2.18.0 and 2.18.1)

The issue is that when I run this script (I have pared it down to the essentials for this case that still causes error)  I often (but not every time!! ) get the error:

  python.exe has stopped working

Note that the script still operates as intended!!!  The script works!! but this error causes the script to pause in the middle requiring a click of the OK button.

Also, at the end of the script I get an error:

  wrapped C/C++ object of type QProgressBar has been deleted See log for
  more details

Log message:

  2016-10-01T12:42:36   2   Uncaught error while executing algorithm
              Traceback (most recent call last):
                File "C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 204, in execute
                  progress.setPercentage(100)
                File "C:/PROGRA~1/QGIS2~1.16/apps/qgis/./python/plugins\processing\gui\MessageBarProgress.py", line 55, in setPercentage
                  self.progress.setValue(i)
              RuntimeError: wrapped C/C++ object of type QProgressBar has been deleted

This error happens every time.

I think that this is happening more when the rasters I am operating on are large
(I am merging 4x 40 Mb rasters in this case)

Here is the excerpt of my script pared down to show the actual activities that still cause consistent failures.

I am not sure if this is a problem with my code or something else!

from qgis.core import *
from qgis.utils import *
from qgis.gui import *
import processing
from PyQt4.QtCore import *
import time

proj = QgsProject.instance()
proj.clear()

workdir = "C:/QGISwork/ving_script/"
vindef = "FR_83510_00001"
date = "20160830"

proj = QgsProject.instance()
proj.clear()

curfield =1

GRN_raster = workdir + vindef + "_"+ str(curfield).zfill(2) +"_" +date +"_"+ "GRN" + ".tif"
RED_raster = workdir + vindef + "_"+ str(curfield).zfill(2) +"_" +date +"_"+ "RED" + ".tif"
RDE_raster = workdir + vindef + "_"+ str(curfield).zfill(2) +"_" +date +"_"+ "RDE" + ".tif"
NIR_raster = workdir + vindef + "_"+ str(curfield).zfill(2) +"_" +date +"_"+ "NIR" + ".tif"

merge_rast = "C:/QGISwork/ving_script/merge_int_"+str(curfield).zfill(2)+".tif"
print "Field number: " + str(curfield).zfill(2)
print "trying merge"
processing.runalg('gdalogr:merge', GRN_raster + ";" + RED_raster + ";" + RDE_raster + ";" + NIR_raster, True, True, 5, merge_rast)

print "Merge complete"


texts: Shapely coordinate sequence to GeoDataFrame

texts:
Using the shapely bounding box function box, and a GeoDataFrameof points, I have extracted a bounding box and created a buffer around it (below gdf is a GeoDataFrameof points):
bbox=shapely.geometry.box(*gdf.geometry.total_bounds)
bbox=bbox.buffer(5)

Now I check and see that bbox is a shapely polygon:
type(bbox)
shapely.geometry.polygon.Polygon

My goal is, from this shapely polygon, to create a 1-row GeoDataFramewith Polygon geometry.  I know that I can get all the coordinates that make up this polygon with:
[*bbox.exterior.coords]

This returns the list of coordinates from a coordinate sequence object which looks like this:
[*bbox.exterior.coords][:5]
[(151.20138500000007, -43.787852999999984),
 (-151.40006999399998, -43.787852999999984),
 (-151.89015569564776, -43.76377663336097),
 (-152.3755216040806, -43.69177940201614),
 (-152.85149338027227, -43.57255467866103)]

However, I cannot seem to figure out how to convert these pairs of coordinates into a GeoDataFrame.  I have tried:
crs = {'init': 'epsg:4326'}
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[*bbox.exterior.coords])

which returns the error:
TypeError: Input must be valid geometry objects: (151.20138500000007, -43.787852999999984)

How, then, does one make these coordinate pair tuples into valid geometry objects for a GeoDataFrame?


texts: PostGIS with clause not using spatial index

texts:
This query uses the spatial index and it's quite fast:

 explain
    select a.the_geom from
     chicago_major_streets a,
     stormevents_locations_2014 b
    where st_intersects(a.the_geom,b.the_geom)

  QUERY PLAN
------------------------------------------------------------------------------
 Nested Loop  (cost=0.28..32436.92 rows=88 width=137)
   ->  Seq Scan on chicago_major_streets a  (cost=0.00..1054.65 rows=16065 width=137)
   ->  Index Scan using stormevents_locations_2014_the_geom_geom_idx on stormevents_locations_2014 b  (cost=0.28..1.93 rows=2 width=32)
         Index Cond: (a.the_geom && the_geom)
         Filter: _st_intersects(a.the_geom, the_geom)
(5 rows)

However if I use a with clause index are ignored:

explain
    with a as (select * from chicago_major_streets),
      b as (select * from stormevents_locations_2014)
    select a.the_geom from
     a,
     b
    where st_intersects(a.the_geom,b.the_geom);

QUERY PLAN
----------------------------------------------------------------------
Nested Loop  (cost=2534.91..232135969.50 rows=283954 width=32)
 Join Filter: ((a.the_geom && b.the_geom) AND _st_intersects(a.the_geom, b.the_geom))
 CTE a
   ->  Seq Scan on chicago_major_streets  (cost=0.00..1054.65 rows=16065 width=409)
 CTE b
   ->  Seq Scan on stormevents_locations_2014  (cost=0.00..1480.26 rows=53026 width=113)
 ->  CTE Scan on a  (cost=0.00..321.30 rows=16065 width=32)
 ->  CTE Scan on b  (cost=0.00..1060.52 rows=53026 width=32)
(8 rows)

Any ideas?


texts: How to see previous years' satellite imagery on GEE?

texts:
I have a set of points for which I want to see what built up or non-built areas were around these points in 2014, 2015, 2017 and 2018. When I lay these points (in yellow) on the default satellite imagery that appears in GEE, I get this:

Now I can zoom in on the map around these points and see actual buildings, roads or fields. Likewise I want to see these points on a similar satellite imagery of 2014, 2015, 2017 and 2018.  I tried using image from Landsat 8 Raw Scenes (Orthorectified) as well as Landsat 8 C1T1 and Real-time data raw scenes, both with their rgb bands. In both of them, we can see the difference in terrains but the details of buildings, roads and fields are lost. I wanted to know if there is any date/year field or some other Google Earth Engine APIs I can use to get the similar map just for the previous years and not the most recent one?


texts: Converting a Multipart Polygon to a Single Line

texts:
I am dealing with a multipart polygon shapefile of Canada, where each province is its own polygon. Simply for cartographic purposes I want to display the provincial boundaries as a dashed line but given that shared boundaries are not actually shared (2 lines overlapping) in a multipart polygon I am unable to use a dashed line as it does not line up between the two shapes and is no longer dashed!

My end goal would be to 'merge' those shared boundaries so they appear as a single line rather than two overlapping lines. I have tried dissolving, multipart to singlepart, topology, polygon to polyline and nothing seems to achieve my goal.

Like I said, this is purely for cartographic purposes so the method does not have to be perfect and a MacGyver workflow would suffice. I am running ArcGIS 10.4.1 for Desktop basic lisence.


texts: How to disable imageSmoothingEnabled for geoImageLayer?

texts:
Considering the code below (which uses openlayers and Viglino's ol-ext), how do I disable imageSmoothingEnabled for the geoImageLayer? I am trying to achieve the same thing as in this question, but for the geoImageLayer, which is part of Viglino's ol-ext.

var map = new ol.Map({
  target: 'map',
  view: new ol.View({
    zoom: 18,
    center: ol.proj.fromLonLat([4.8728,52.3316])
  }),
  layers: [new ol.layer.Tile({
    source: new ol.source.OSM()
  })]
})

var geoImageLayer = new ol.layer.Image({
  name: "Georef",
  opacity: 0.5,
  source: new ol.source.GeoImage({
    url: 'image.jpg',
    imageCenter: ol.proj.fromLonLat([4.8728,52.3316]),
    projection: 'EPSG:3857',
    imageScale: [1, 1]
  })
})
map.addLayer(geoimg);

The internet says that this should work:

map.on('precompose', function(evt) {
  evt.context.imageSmoothingEnabled = false;
});

And it does work for the background OSM layer, but it does nothing for the geoImageLayer

I've also tried,

geoImageLayer.on('change', function(evt) {
  evt.target.values_.source.canvas_.canvas_.getContext('2d').imageSmoothingEnabled = false;
});

but it didn't do the trick. Probably because I'm modifying a copy of the actual layer, but I don't have any other ideas.

Here is a jsfiddle with the above code.


texts: Geometry not displaying correctly in QGIS 2.12

texts:
I loaded a few Shapefiles from different sources into a PostgreSQL database, when I visualize them in QQGIS some of them are not shown.

This how i added the geometry column in my tables:

SELECT AddGeometryColumn ('public','table_geom','geom',4326,'MULTIPOLYGON',2, true);

And this is the SQL script that I used to change the srid of the Shapefiles I uploaded to the table:

SELECT UpdateGeometrySRID('public', 'shapefile_table','geom',4326);

I added the geometry column from the table created from the Shapefile, there it wasn't any trouble with this.

This is how i defined my table:

  CREATE TABLE table_geom
    (
      id serial NOT NULL,
      geom geometry(MultiPolygon,4326),
      CONSTRAINT pk_table_geom PRIMARY KEY (id_ageb),
      CONSTRAINT enforce_dims_geom CHECK (st_ndims(geom) = 2),
      CONSTRAINT enforce_geotype_geom CHECK (geometrytype(geom) = 'MULTIPOLYGON'::text OR geom IS NULL),
      CONSTRAINT enforce_srid_geom CHECK (st_srid(geom) = 4326)
    )
    WITH (
      OIDS=FALSE
    );
    ALTER TABLE table_geom
      OWNER TO postgres;

The thing is, i can visualize correctly the source Shapefile in QGIS  but not from the database.
What am I missing? I already checked the srid is the same in all tables but that seems that is not the problem.
Any hints?

Quick edit: i think is a bug(?), probably because the table data is too big, now some of the layers are being shown from one table, but not all.

Edit 1: I see now what is the problem, the geometry from the table is being drawn but not correctly,as i mentioned before if i load the source shapefile in qgis the area is shown where it should be, but if i load from the database this area is drawn elsewhere, but why? i changed the srid in all loaded shapefiles


texts: KnownCoordinateSystems.Geographic.Europe.OSGB1936 showing as WGS84 in Visual Studio

texts:
I'm trying to get my head around DotSpatial.
I'm specifically trying to convert the geographic points of a feature from OSGB1936 (27700) to WGS84.
This is the line of C# that I was hoping would do the conversion:
Reproject.ReprojectPoints(polygonPoints, z, KnownCoordinateSystems.Geographic.Europe.OSGB1936,
KnownCoordinateSystems.Geographic.World.WGS1984, 0, feature.Coordinates.Count);

The coordinates after conversion are incorrect (all are the value 8).
In Visual Studio if I hover over KnownCoordinateSystems.Geographic.Europe.OSGB1936 I get "{WGS84}" returned from Visual Studio suggesting it has evaluated that expression to that coordinate system. DotSpatial is recognising KnownCoordinateSystems.Geographic.Europe.OSGB1936, i.e. it is not highlighting as an error and stopping execution, but it seems to be ignoring it and returning WSG84 instead.
I've tried installing other NuGet DotSpatial packages like DotSpatial.Projections and have searched the web but cannot find out why. I feel like I'm missing some definitions somewhere.
Can anyone help me out?


texts: Storing gdal2tiles output directory in MongoDB, add as XYZ layer in OpenLayers

texts:
I am regularly running gdal2tiles.py which generates millions of XYZ image tiles for me. Right now I either store these tiles as flat files on a Windows Server file server or I upload a directory to Google cloud storage if needed.

I use a virtual directory or access the tiles directly from cloud storage to display XYZ tile layers in OpenLayers.

I would like to test storing the tile directories in MongoDB. Has anyone done a script to store the gdal2tiles.py output directly in MongoDB?

Or has anyone a script to upload an existing file directory into MongoDB?

Does anyone use MongoDB to store image tiles?

I would like to compare the performance with my traditional file server and Google cloud storage.


texts: Read Postgis geom available as string variable with R

texts:
I have a Postgis Database from which I query data. I know that I can use the rgdal::readOGR() function with the PG driver to get geoms from the database into R.

However, I am also querying lots of non-spatial data from the same database using DBI. From a programming point of view it seems not clean to have two different methods to connect to the same database in the same script.

What I do now is to query the column with SELECT ST_AsBinary(the_geom) ... using DBI and then read the result with wkb::readWKB(). This works (only uses DBI for connection) but is quite slow.

Is there any way to read a geom string as returned from Postgis with rgdal when the string is already available as variable?


texts:
How do I convert from an ESRI.ArcGIS.Geometry.IPoint to an ESRI.ArcGIS.Client.Graphic point object?

In other words, I'm wanting to convert from an IGeometry to a runtime geometry, I guess.


texts: Correct Buffers crossing dateline in sf package in R

texts:
Using the global coral reef distribution data provided by UNEP https://data.unep-wcmc.org/datasets/1

I have been attempting to buffer around the points and polygons of coral reefs. The point buffers (green) seem to be fine, however the polygon buffers (pink) are seemingly spread across the world.

I was wondering if there is a technique to clean up the buffers that spread across the world? Or is there geometries with the data that is causing this to happen? I am not sure I am asking the right question here.

Here is my code for the multi polygon data, which I replicated for the point data:

##projection,

c_py <- st_transform(c_py, crs = "+proj=laea +lat_0=45.5 +lon_0=-114.125 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")

##buffers 100km

buffer_100_py <- st_buffer(c_py, dist = 100000)

After which I printed a map using tmap package, where you can see the buffers in pink spread across the whole map.

Here is the buffers after st_union, if this may reveal better what behaviour is going on here?

I am new to spatial data, and was wondering if there was something simple I was missing?

Update:

In an attempt to correct for the buffers crossing the dateline I applied this code to use the st_wrap_dateline function.

buffer_100_py2 <- buffer_100_py %>% st_wrap_dateline(options = c("WRAPDATELINE=YES",  "DATELINEOFFSET=180"), quiet = TRUE)

Note here is just polygon buffer, point buffer excluded

So can see the function hasn't worked as I expected.


texts: Print composer forgets layers and puts them in wrong order

texts:
I have a project that mainly consist of PostGIS tables, but also some shapefiles. All in all I have some 50 layers in use in the project.

The last days, I have had problems with the print composer. Two polygon layers are shown behind point and line layers even if layer ordering is enabled (the same problem in exported maps), they are ordered correctly in QGIS desktop, and shown in the right order in QGIS desktop. Also, at least one line layer does not show in the print composer and is not exported when creating an exported map.

I have tried removing and importing PostGIS layers, same thing happens again. I have also tried removing layer styles and setting them up again without any extra settings.

This is QGIS 2.18. I cannot move to QGIS 3 because of a geomnearest function that does not work with QGIS 3, otherwise that would be my next thing to try.

Are there any ways to check the coherency of the project file? Is there anything else that comes to mind that could cause this kind of trouble?

This has something to do with visibility presets. Even if I have the correct layers saved as a visibility preset, and I tell the composer map to follow that visibility preset, some layers are not shown in the composer (but shown alright in QGIS Desktop when the visibility preset is chosen there). When i deselect the use of visibility preset in the composer, the layers show OK. I have tried to delete presets and re-make them with new names, does not help.


texts: I have generated a delaunay polygon file that represents terrain data. Is there a way I can produce contours from the polygon file? I know I can convert the polygon to an Esri TIN, and then produce other products, but is there a way to do it directly from the polygon?


texts: What is the way to calculate center possitions from WKT format POLYGON?


texts: Capturing Output Messages from GRASS when Running Python Scripts

texts:
Using documentation here:
https://grasswiki.osgeo.org/wiki/Working_with_GRASS_without_starting_it_explicitly#Python:_GRASS_GIS_7

I've figured out how to access GRASS modules in Python without explicityly starting GRASS; I set up the environment and run my script in IDLE and all goes well. I'm using GRASS 7 on Linux with Python 2.7.

What I haven't been able to figure out is how to capture the output that GRASS would generate and print to the command console. I would like to print that information to the screen so that I have an indication that things worked. For example, running this command from either the GRASS shell or gui:

gscript.run_command('v.in.ascii', input=pntfile, output='raster_pnts',
                  overwrite=True)

Would print this to the console and do the operation:

  (Fri Jan 15 16:17:35 2016)
  v.in.ascii --overwrite input=/grass_work/output/points_out.txt
  output=raster_pnts WARNING: Vector map  already exists
  and will be overwritten Scanning input for column types... Number of
  columns: 3 Number of rows: 93 Importing points... Populating table...
  Building topology for vector map ...
  Registering primitives... 93 primitives registered 93 vertices
  registered Building areas... 0 areas built 0 isles built Attaching
  islands... Attaching centroids... Number of nodes: 0 Number of
  primitives: 93 Number of points: 93 Number of lines: 0 Number of
  boundaries: 0 Number of centroids: 0 Number of areas: 0 Number of
  isles: 0 (Fri Jan 15 16:17:36 2016) Command finished (0 sec)

When I execute my Python scripts I also want this to print to the screen, but I can't figure out how to do it. I know how to do this with certain commands where strings are returned as output, like:

region = gscript.parse_command('g.region', flags='p')
print region

Would parse and print the region info for me. I know this is different, as I'm trying to capture error messages or progress and not the output itself. I have experimented with some of the examples here:

https://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library

And have tried using start_command instead of run_command to try and capture and print output messages, but I can't get it right (have experimented and searched for several hours).

(I also tried using subprocess.check_output - I've successfully used that for passing input from Python to the GDAL/OGR command line tools and getting the output messages back, but can't get it working here).


texts: How to create a human migration map with QGIS

texts: I want to create a Human migration map looks like this http://postimg.org/image/an217oxh9/ or http://postimg.org/image/6ryn54ebh/ using QGIS, but don't know how to that because arrows available in QGIS are straight arrows while I need curved arrows. What plugin should I use to have such curved arrows.


texts: How to install Colour Scale Bar for One-Band Rasters

texts: I am new to QGIS. I want to generate a Colormap of a picture. I searched and found one plugin-Colour Scale Bar for One-Band Rasters-can do this. But I didn't find it the the plugins. Is anyone know why?


texts:
I created an atlas in composer with my results from site visit forms.

I would like to display a checkbox with ticked or unticked state depending on my field value.

I am thinking as a workaround to use two svg files which are displayed or not depending the field value, but it is a bit complicated when multiplicating the fields.

Is there an easier way to do so?


texts:
I am trying to join a shp and dbf file using the PyQGIS API.
The code I'll use is below:

from PyQt4.QtCore import *
vector1 = QgsVectorLayer("path to dbf file", "test1", "ogr")
id1 = vector1.id()
print id1
vector2 = QgsVectorLayer("path to shapefile", "test2", "ogr")
id2 = vector2.id()
print id2
QgsMapLayerRegistry.instance().addMapLayer(vector1)
QgsMapLayerRegistry.instance().addMapLayer(vector2)
info = QgsVectorJoinInfo()
info.joinFieldIndex = QgsMapLayerRegistry.instance().mapLayer(id1).fieldNameIndex("Id1Column0") # Returns 0 as expected
info.joinFieldName = QString(u"Id1Column0")
info.joinLayerId = QgsMapLayerRegistry.instance().mapLayer(id1).id()
info.targetFieldIndex = QgsMapLayerRegistry.instance().mapLayer(id2).fieldNameIndex("Id2Column0")  # Returns 0 as expected
info.targetFieldName = QString(u"Id2Column0")
info.memoryCache = True
print info.joinFieldIndex #All the following show the correct values
print info.joinFieldName
print info.joinLayerId
print info.targetFieldIndex
print info.targetFieldName
print info.memoryCache
QgsMapLayerRegistry.instance().mapLayer(id2).addJoin(info)

Querying for the respective field names using QgsVectorLayer.attributeDisplayName(0) yields the respective column names as expected.

Yet when I look in the GUI properties, while the join layer is set correctly, the join field is something different and the target field is empty.

Could anyone shed some light on what goes wrong?


texts: Binding to Selected Feature ESRI/WPF Runtime/MVVM

texts:
Given a ArcGISLocalFeatureLayer as a layer on the WPF runtime map object, how can I bind a ViewModel to the selected feature or graphic?

<esri:ArcGISLocalFeatureLayer x:Name="MySampleLayer" ID="MyFeatures"
Path="C:/..../Samples.mpk" LayerName="MyFeatures" OutFields="*" DisableClientCaching="True"
Editable="True" AutoSave="False"  ValidateEdits="True"
Mode="OnDemand"  SelectionColor="#FFFFFF00" />

I realize there is already a FeatureDataGrid and FeatureDataForm in the runtime, but I'm trying to achieve a little different user experience. For simplicity, suppose I just have some controls on the same window as the map:

<Grid >
    <StackPanel>
        <StackPanel.DataContext>
            <local:FeaturePropertiesViewModel/>
        </StackPanel.DataContext>
        <TextBox Text="{Binding FeatureName}" />
        <TextBox Text="{Binding FeatureInfo}" />
    </StackPanel>
</Grid>

and a simple ViewModel ...

class FeaturePropertiesViewModel
{
    public string FeatureName { get; set; }
    public string FeatureInfo { get; set; }

    public Graphic SelectedFeature { get; set; }     // ??? use this somehow?
    public GraphicCollection AllGraphics { get; set; }  // ??? use this somehow?
}

How do I get those text fields updated to when a new feature is selected? Ideally, I'd like two-way binding so that changing the text updates the underlying feature class.

Update: Based on Darlene's suggestion, I started looking into other ways to get the information include the mouse interactions - since that is how the feature is being selected.

I can capture/handle the MouseLeftButtonDown event and do everything I need to in the code-behind. So I started working of that to refactor to MVVM. Using MvvmLight, I can map that event to a command like this:

<esri:ArcGISLocalFeatureLayer ... />

            <i:Interaction.Triggers>
                <i:EventTrigger EventName="MouseLeftButtonDown">
                    <cmd:EventToCommand Command="{Binding Path=MyLeftMouseCommand,
                     Source={StaticResource MyVieModel}}"
                     CommandParameter="{Binding  RelativeSource={RelativeSource AncestorType=esriLocal:ArcGISLocalFeatureLayer}}"
                     PassEventArgsToCommand="False" />
                </i:EventTrigger>
            </i:Interaction.Triggers>
</esriLocal:ArcGISLocalFeatureLayer>

However, I need both the mouse event args (GraphicMouseButtonEventArgs) and the selected graphic. And I'm having trouble getting both passed to a command. I tried RelayCommand like this ...

public RelayCommand<ArcGISLocalFeatureLayer> MyLeftMouseCommand{ get; protected set; }

In constructor ...

    MyLeftMouseCommand= new RelayCommand<ArcGISLocalFeatureLayer>(p =>
    {
        Console.Writeline("It works!");
    });

And that will get me the CommandParameter value if I set PassEventArgsToCommand="False", or I can get GraphicMouseButtonEventArgs if I set PassEventArgsToCommand="True", but I can't get them at the same time!

Also using ReactiveUI if there is anything in there that might be helpful.


texts:
I followed the instructions here:

How to install QGIS on Centos 6?

If I run:

QgsVectorLayer * vpoly = new QgsVectorLayer("Polygon", "pointbuffer", "memory");
QgsVectorDataProvider * provider = vpoly->dataProvider();

Then provider is a NULL pointer

Also, if I do:

QString myPluginsDir = "/usr/lib64/qgis";
QgsProviderRegistry * preg = QgsProviderRegistry::instance(myPluginsDir);
QString pluglist=preg->pluginList();
printf("plugins: %s\n",pluglist.toStdString().c_str());

Then it prints:

plugins: No data provider plugins are available. No vector layers can be loaded

Also, when I try to load a QgsVectorLayer with the following:

QgsVectorLayer * mypLayer = new QgsVectorLayer(myLayerPath, myLayerBaseName, myProviderName);
if (mypLayer->isValid()){
    qDebug("Layer is valid");
  } else {
    qDebug("Layer is NOT valid");
    return;
  }

Then it says that the layer is NOT valid.

This all works just fine in Ubuntu, but I can't get it to work in centos.  I think maybe it's missing all the plugins.  What am I missing?

Please help.