texts: ESRI Silverlight Cartography - Looking at whats not there

texts:
I have line data that has been symbolozied in ways that I don't have time to understand and would like to make a map that tells me =).

My general question is how to display the lack of data results for queries.  I have a bunch of layer files that were built with a bunch of def queries in ArcMap.  I will iterate over the line data and for each unique line check for results of applying a def query.  If results are returned, a new layer file is created representing that query.

I want to highlight the queries that didn't return results in a silverlight app.  So I have all the queries that returned zero records in a database.

I'll be making this map with silverlight but any solid cartography concepts for getting at what's not there are welcome.


texts:
I cannot see how Union is being used; why in this instance and whether it produces the desired results.
I am posting the procedure for some clarity. I don't understand the selection its asking for either; could someone maybe walk me through that as well?
I know I'm blanking on something simple but none of the bosses even understand GIS
Percent Impervious of Managed Drainage Areas
Calculate the % of the managed area covered by impervious surface for each project.
Required Data:
»»
GSI_Projects
»»
GSI_DrainageAreas (only those associated with GSI footprint locations)
»»
Impervious Area
PROCEDURE:

Add 2 fields “PercDAImp” and “ImpDAsq” to the GSI_Projects file

**2. Create a union between the GSI_DrainageAreas and the Impervious Surface layer.

Add Field “Impervious” to the unioned file

Select all rows that intersected both the drainage areas and impervious

Calculate Impervious using Field Calculator as “Yes”

Switch selection and calculate Impervious as “NO” – clear all selections**

Dissolve unioned file by GSI_ID and Impervious

Add 2 fields – ImpDAsf and PervDAsf

Select all Impervious = Yes and use the Calculate Geometry to calculate the square footage into ImpDAsf

Switch selection and use the Calculate Geometry to calculate the square footage in the PervDAsf. Clear all selections

Dissolve by GSI_ID summing fields ImpDAsq and PervDAsf and title
GSI_DrainageAreas_ImperviousTotal

Add a field PercImperv

Calculate field as [ImpDAsq]/([ImpDAsq]+[PervDAsq])

Join the file GSI_DrainageAreas_ImperviousTotal with GSI_Projects by GSI_ID

Use the field calculator to transfer the PercImperv as well as the ImpDAsq fields into GSI_Projects.

The proportion of impervious area managed will serve as the score for this criteria


texts: How do I search datasets through the using search using our Enterprise?

texts:
So when you have layer that is missing or has a exclamation point next to the layer.
Now, I am not sure where to find the dataset because our enterprise system in Citrix has thousands of dataset all over.
When I use that search toolbox to find the layer that is missing and it came up nothing or doesnt work.
It would say "Warning  Enterprise search services are not found.
Use Search Options dialog in the search toolbar to register them. I have never use that system in Citrix.
Can you please advise me?
We have C hard drive on our computer but the T drive is what we use to find or save our datasets on our network.

I am using ArcGIS 10.2.2 and using Citrix.


texts: Lat/Long attribute to Trimble Pathfinder Data Dictionary

texts: Is it possible to add beginning and end point lat/long attribute to Trimble Pathfinder Data Dictionary before collecting data in the field?


texts: Running python script with HTML and Javascript on Map Tips

texts:
I have a problem of running python script by pressing a HTML button on my Map Tips.

Right now, this is what I've written for my Map Tips in HTML format.

<script>
function option1(){
    [% testingFunc('test1') %]
}
function option2(){
    [% testingFunc('test2') %]
}

<button onclick="option1()">Option #1</button>
<button onclick="option2()">Option #2</button>

My aim is to allow data input by pressing button that shows up in the Map Tips (http://i.imgur.com/lrzMT6O.png). Each button will run a function named "testingFunc" with different input.
Currently, for testing purpose, testingFunc only prints out the input of the function to the console and does nothing else.

Unfortunately, what happens when I hover my mouse over the points (to display the Map Tips,) it ran both scripts immediately and when I pressed the buttons, they don't do anything at all.

I'm guessing this is because the QGIS expression takes precedence while executing that HTML code, but that makes me wonder if what I wanted to do is even possible to begin with?


texts: Combining the output of ST_Collect/ST_Union to create linestring from multilinestring?

texts:
I am using pgr_dijkstra() to compute the shortest path and then using ST_Collect (or ST_Union) to get the shortest path geometry. However, the collection returned from the ST_Collect is a multi-linestring and not amenable to my further analysis which includes using functions like ST_LineLocatePoint etc. which need a linestring as an argument.

I am having to use ST_Collect() instead of ST_MakeLine() directly on the pgr_dijkstra output as ST_MakeLine ends up connecting the sub-geometries out of order resulting in a garbled line (described here).

Function code:

CREATE OR REPLACE FUNCTION public.sp_od(
    orig integer,
    dest integer)
    RETURNS TABLE(shortest_path geometry)
LANGUAGE 'sql'

AS $BODY$

SELECT st_makeline(geom) as shortest_path
FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom, true) as cost FROM public."WA_roads"',
(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(ocentx as double precision), CAST(ocenty as double precision)), 4326) from all_trips_non_zero where origin = orig LIMIT 1) ASC
LIMIT 1),
(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(dcentx as double precision), CAST(dcenty as double precision)), 4326) from all_trips_non_zero where destination = dest LIMIT 1) ASC
LIMIT 1), directed := false
) as pt
JOIN public."WA_roads" rd ON pt.edge = rd.id;

$BODY$;

I tried using ST_SnapToGrid() and ST_LineMerge(), still the resulting output from ST_Collect is a multilinestring.

I seek suggestions to either get the output from pgr_dijkstra() in the right order or for converting the ST_Collect() multilinestring to linestring.


texts:
You know you can use CQL to create GeoTools Filter used when filtering FeatureCollection. GeoTools 16 user guide has made a quick tutorial and examples for CQL strings like the following ones:

ATTR1 >= 5
NCLUDE
EXCLUDE
ATTR1 < (1 + ((2 / 3) * 4))
(ATTR1 + ATTR2) > 0.07
ATTR1 < abs(ATTR2)
ATTR1 < 10 AND ATTR2 < 2 OR ATTR3 > 10
ATTR1 LIKE 'abc%'
ATTR1 NOT LIKE 'abc%'
ATTR1 BETWEEN 10 AND 20
ATTR1 IS NULL
ATTR1 IS NOT NULL
ATTR1 DOES-NOT-EXIST     // check if something does not exist

string functions:

strToLowerCase(ATTR1) like '%m%'

collections:

ATTR1 IN ('New York', 'California', 'Montana', 'Texas')
IN ('states.1', 'states.12')       // with specific feature ids

date time:

ATTR1 TEQUALS 2006-11-30T01:30:00Z                       // equal to a date
ATTR1 BEFORE 2006-11-30T01:30:00Z                        // before a date
ATTR1 BEFORE 2006-11-30T01:30:00Z/2006-12-31T01:30:00Z   // before a period
ATTR1 AFTER 2006-11-30T01:30:00Z                         // after a date
ATTR1 AFTER 2006-11-30T01:30:00+03:00                    // after a date using time zone GMT+3
ATTR1 AFTER 2006-11-30T01:30:00Z/2006-12-31T01:30:00Z    // after a period

geometry:

CONTAINS(ATTR1, POINT(1 2))
BBOX(ATTR1, 10,20,30,40)
DWITHIN(ATTR1, POINT(1 2), 10, kilometers)
CROSS(ATTR1, LINESTRING(1 2, 10 15))
INTERSECT(ATTR1, GEOMETRYCOLLECTION (POINT (10 10),POINT (30 30),LINESTRING (15 15, 20 20)) )
CROSSES(ATTR1, LINESTRING(1 2, 10 15))
INTERSECTS(ATTR1, GEOMETRYCOLLECTION (POINT (10 10),POINT (30 30),LINESTRING (15 15, 20 20)) )

full list of geometric predicates (from geoserver guide ):

EQUALS, DISJOINT, INTERSECTS, TOUCHES, CROSSES, WITHIN, CONTAINS,
OVERLAPS, RELATE, DWITHIN, BEYOND

examples are usefull but I couldn't find any full tutorial about syntax, structure, functions, math functions, etc. Can anyone introduce a full guide for CQL.

udig tutorial:

uding guide for CQL


texts: Interpolate nodata in rasters using Python

texts:
I have a raster of the shape (1000,1000) and some areas having no data values. I would like to fill the data gaps by interpolating or tinning (does not matter) over the surrounding areas, however I fail to do that using Python.

I have searched tried some procedures already discussed at stackexchange, but failed to succeed:

Attempt 1: using gdal.FillNodata

import gdal, gdalconst
file = "C:\\Foo\\bar.tif"

raster = gdal.Open(file, gdalconst.GA_Update)
fill = raster.GetRasterBand(1)
f = gdal.FillNodata(targetBand = fill, maskBand = None, maxSearchDist = 50, smoothingIterations = 0)

fill= fi.ReadAsArray()
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create("D:\\Foo\\out.tif", cols, rows, 1, gdal.GDT_Float32)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(fir)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()

### does not work - it takes ages, even with a smaller example. And if it works, there are no changes compared to the input file

Attempt 2: Using rasterio.fill -- However, this returns an AttributeError.

Attempt 3: Using Scipy

import tifffile
from scipy.interpolate import griddata

raster = tifffile.imread("D:\\Foo\\bar.tif")

grid_x, grid_y = np.mgrid[0:1000, 0:1000]

nans = np.argwhere(np.isnan(raster))
coordnotnan = np.argwhere(~np.isnan(raster))
valuenotnan = excr[coordnotnan]
vanona = valuenotnan[:,0,0]

grid_z1 = griddata(coordnotnan,vanona,nans, method=nearest) #for sake of testing

grT = grid_z1.T
rR= list(raster.ravel())
new = []
i = 0
for na in rR:
    if np.isnan(na) == True:
        new.append(grT[i])
        i += 1
    else:
        new.append(na)

arrayed = np.array(new, dtype=np.float32).reshape((1000,1000))
#...
## Besides being inefficient (loops), it does not seem to take the values of the neighboring pixels

Do you spot the mistakes I've made? Is there another way to get this job done?


texts: Raster Dataset Misaligning in ArcGIS Desktop?

texts:
We have an ArcSDE Raster Dataset of aerial imagery covering a County. A little while back the dataset somehow got tweaked to where it doesn't align correctly anymore.  It is in the NAD83 StatePlane coordinate system.

We could georeference the image but we're not sure how it got misaligned in the first place, therefore, perhaps the solution is something other than georeferencing?

I am using ArcMap 10.6.


texts:
Is there a tool that does the opposite of a dissolve?

I would like to break a shapefile into multiple files based on an attribute.


texts: Adding parameter to all requests made by ArcGIS iOS SDK?

texts:
I am developing an app that uses ArcGIS iOS SDK.

My requirement is that I want to pass a session token (Received from my backend along with Login response) to each and every service call the SDK makes. As we know the SDK exposed API & usually accepts URL as a parameter. If i supply the URL of the map service, it would fetch the map info like images, etc & loads it on my view. This token is used by my backend team to authorize the user & respond the relevant map information from GIS server.

I know about AGSCredential object that is used to access secure GIS servers. My requirement is not only connecting to a secured GIS server but also get data based on current user. My app has 5 kinds of users who access the same GIS but should not be able to access complete data. To achieve this the Backend team has introduced a proxy layer, where the iOS SDK should talk to it. The proxy layer responds with appropriate data based on token that the client sends along with every request. The SDK requests should have a header with the token.

I wanted to know if this possible and, if yes, how can I accomplish it?

The web developers who use the ARCGIS SDK for .Net managed it by writing a proxy between the SDK and service layer on client. They are suggesting me to also do a similar thing but I doubt if that is possible in iOS.

Perhaps someone came across a similar scenario where you need to proxy the service calls made by a library integrated in an iOS application. I am sure we can add a parameter or header if the service request is written by me, but I want to proxy the calls happening from my SDK (ArcGIS iOS).


texts: QGIS r.water.outlet multiple points processing modeler

texts:
I am trying to create a process with the process modeler. I want to run r.water.outlet by importing multiple points (a shapefile of points).
I searched through the internet and StackExchange. Could you help me automate my process?


texts: Identify corrupt shapefiles or rasters via Python or ArcObjects

texts:
I'm interested in a way to identify rasters (e.g., ESRI GRID format) or shapefiles that are corrupt that would prevent them from being used in an analysis tool. I have a large number in structured folder tree and would like to be able to script something that will walk the file tree and be able to ID rasters that, for example if you tried to view it in ArcCatalog it would not allow it to be displayed and accordingly would make any geoprocessing tool fail. Same deal for shapefiles or feature classes.

My preference would be somehow to do this in Python (and I wouldn't be opposed to a non-ESRI python module or other) but I'm open for anything. I guess you could hack something together by trying to validate geometry or capitalizing on some aspect of the file design, but what would make the most sense? Or is there another way to do it that is specifically designed for such a purpose?

My goal is to just be able to run a simple Quality Control mechanism on some large week-long preprocessing files to make sure they're all sound before I start doing the actual batch processing.

thanks,
Tom


texts: How to call GEOS error handlers provided during initialization (initGEOS)?

texts:
From what I learned, when using GEOS in C/C++, it should be initialized first using initGEOS :

extern void GEOS_DLL initGEOS(GEOSMessageHandler notice_function,
         GEOSMessageHandler error_function);

, which takes two error handling functions notice_function, error_function.

My question is:

How can these two handler functions be accessed and called in C? (e.g. when one wants to print some notice messages while debugging an existing GEOS function)

I can't seem to find the documentation about this.


texts: Is QGIS unable to modify layers from a .gdb?

texts:
There is a .gdb file with several census tables I downloaded from TIGER. I am trying to join one of the tables to a shapefile but am not able to as the join field for the table contains excess characters. When I go into the field calculator to make a new attribute field without those characters, I found I was not able to and was greeted with: "This layer does not support adding new provider fields."

Is it possible to somehow edit this attribute field in QGIS or will I need to use ArcGIS for it?


texts:
I've been looking at some various ways to calculate geodesics. One particularly good site collecting load of helpful information, is this one. There you can also find a link to a map that shows a circle whose geodesic radius is accurately calculated by the WGS84 ellipsoidal approximation. In the resulting map, there are some additional info presented (as red arcs), called envelope. I'm struggling to understand what these are and what exactly they represent. The only partially useful resource I can find is from this website, but looking at the faint lines drawn on 2D image of a stone, does not enlighten anyone, especially as you cannot rotate the images, to see where those lines goes, and intersects.

OR from another calculation:

From the website, the description only say:

  The geodesic envelopes as red curves; all the geodesics emanating from
  lat1, lon1 are tangent to the envelopes (providing they are extended
  far enough). The number of solutions to the inverse problem changes
  depending on whether lat2, lon2 lies inside the envelopes. For
  example, there are four (resp. two) approximately hemispheroidal
  geodesics if this point lies inside (resp. outside) the inner envelope
  (only one of which is a shortest path).

Q: How can I better understand what those red arcs represent?

Why are they useful on this picture?
How are they generated?

I'm more of a visual person, seeking a visual explanation, rather than a purely verbal description.


texts:
I try to set and go to search by help GIS but I couldn't find an answer

How to show coordinate on the polygon but I don't show it on dataframe


texts:
I already have a method to generate greyscale bitmap heightmaps from Google Earth using a long and involved manual process (see here).

Is there a scriptable way of producing a heightmap to a set scale from freely available data, given a set of coordinates?


texts:
I am new to PostGIS.

I have a table of buildings each having a value for number of people living in it. I would like to calculate the total value (sum of number of people) in all buildings which fall into a defined circular buffer.

I have written this code for it which uses a loop, however it gives an error which is followed.

I run the following query in PgAdminIII.

query:

BEGIN
FOR building_pop IN SELECT total_pop FROM Buildings WHERE ST_DWithin(ST_Transform(geom, 26918), (SELECT ST_Transform(ST_GeomFromText('POINT(1116858.6062789 7086290.6680572)', 900913), 26918)), 1000)
    LOOP
    totalpop = totalpop+building_pop;
    NEXT building_pop;
    END LOOP;
    UPDATE public.h_grid200m SET grid_pop = totalpop WHERE gid =1
END

I receive the following error:

  ERROR:  syntax error at or near "FOR" LINE 2: FOR building_pop IN
  SELECT total_pop FROM h_buildings WHERE ...


texts: Collecting and editing local data on a tablet?

texts:
I need to view, collect, and edit data that has to be stored locally on an iPad, though we may be able to get an Android tablet if that's the OS with a solution right now. This data cannot go through the cloud or be hosted on a server.

I've found options like CarryMap that allow you to view data on tablets but nothing that allows you to edit locally. Collector is the perfect app for what we need to do other than accessing local data. Saving it as a tile package would work in some instances but there are times we need more functionality.

Edit: we have data that cannot be stored in the cloud or hosted on a server. It needs to be copied directly to the tablet. Collector can only do that with a tile package basemap but we need more functionality than a static layer.


texts: GeoDjango order_by('distance') not working as expected?

texts:
Using GeoDjango, I am trying to sort a list of users based on the distance from a point. It does generate a list, but its not ordered by distance, though I specify that in my query.
This is the query I am using.
origin = GEOSGeometry('SRID=4326;'+ 'POINT('+ str(buyer.buyer_current_location.x)+' '+str(buyer.buyer_current_location.y)+')')
close_sellers = Seller.objects.exclude(profile__user = user).filter(seller_current_location__distance_lte=(origin, D(mi=distance_mi)),profile__user__is_active = True).distance(origin, field_name='seller_current_location').order_by('distance')[:20]

and these are the models that I have
class Buyer(models.Model):
    # This field is required.
    profile = models.ForeignKey(UserProfile,unique=True)
    # Other fields here
    screen_name = models.CharField(max_length = 200,default='buyer')
    buyer_current_location = models.PointField(null = True, srid=4326)
    objects = models.GeoManager()

class Seller(models.Model):
    # This field is required.
    profile = models.ForeignKey(UserProfile,unique=True)
    # Other fields here
    seller_current_location = models.PointField(null = True, srid=4326)
    objects = models.GeoManager()

Does anyone have any clue regarding what might be going wrong?
Edit
The list seems sorted, but only messes up on for sellers in a certain part of town (its totally crazy!)


texts: Rasterio + QGIS - rasterio.drivers() -> TypeError: 'module' object is not callable

texts:
I have been trying to get rasterio working with my QGIS setup so that I can take advantage of the fast raster manipulations.

After a day of installing / uninstalling and then piece by piece searching and getting components installed one by one where needed (From upgrading GDAL, numpy, and a whole slew of other modules which were fighting each other) I manged to get to the point where I could import rasterio without any errors.

However now when I run, the first call in the example causes a Traceback:

with rasterio.drivers(CPL_DEBUG=True):

creates error =>

TypeError: 'module' object is not callable

I have noticed that other rasterio objects seem fine, but all functions that are in .pyc files have this same issue.

Any help?  - I will continue in the morning, but I was hoping some elves would appear and help solve this one tonight! :-)


texts: Сalculate polygon centroids with 'point on surface' function

texts:
I want to calculate centroid coordinates X_center and Y_center for every polygon in layer so that in case polygon has not simple shape (e.g. with holes or like in the picture below) centroid will be in the polygon anyway:

I found 'point on surface' function in QGIS, that returns a point guaranteed to lie on the surface of a geometry. BUT When trying to calculate coordinates with the function, dump happens, although qgis shows output preview correctly -

I found similar question but the following answer implicates new point layer creation, that is useless for me - I want to get centroid coordinates in the same layer.

Interesting, when creating virtual field - it calculates, I can even pick the polygon and see coordinates, but when opening attribute table, qgis crushes.

What is wrong?


texts: Select the lines that form a ring in postgis

texts:
I need to select all the ids of the linestrings that forms a ring, grouped by ring.

For example, the 6 linestrings below

CREATE TABLE lines(id INTEGER);
SELECT AddGeometryColumn('lines', 'geom', 0, 'LINESTRING', 2);

INSERT INTO lines(id, geom) VALUES
    (1, ST_GeomFromText('LINESTRING(0 0, 1 0, 1 1)')),
    (2, ST_GeomFromText('LINESTRING(1 1, 0 1, 0 0)')),
    (3, ST_GeomFromText('LINESTRING(0 2, 1 2, 1 3)')),
    (4, ST_GeomFromText('LINESTRING(1 3, 0 3, 0 2)')),
    (5, ST_GeomFromText('LINESTRING(2 0, 3 0, 3 1)')),
    (6, ST_GeomFromText('LINESTRING(0 0, 1 0)'));

forms 2 ring: linestrings 1 and 2 form a ring (let's call it "ring 1"), 3 and 4 form "ring2".

I need something like this:

 line_id |  ring_id
 --------+--------------
     1   |     1
     2   |     1
     3   |     2
     4   |     2


texts:
I have a plugin that creates a memory layer and displays it on the canvas.  If I run the plugin a second time after removing the memory layer from the layer list panel without closing QGIS the memory layer is produced along with the previously created layer.  It seems as though the memory layer name cannot be reused until QGIS is closed. Is there a way to clear this layer from the memory to rerun the plugin?

###Create new layer
vecLineURI = "LineString?crs=epsg:4326&field=id:integer"
vecLine = QgsVectorLayer(vecLineURI,"point2line","memory")
pr = vecLine.dataProvider()

###Start editing
vecLine.startEditing()

###Add Feature
feat = QgsFeature()
feat.setGeometry(QgsGeometry.fromWkt(lneWKT))
feat.setAttributes([1])
pr.addFeatures( [feat] )

####Update Extents
vecLine.updateExtents()

###Commit Changes
vecLine.commitChanges()

###Load layer
QgsMapLayerRegistry.instance().addMapLayer(vecLine)

I was removing the layer at the start of the plugin open...

layerMap = QgsMapLayerRegistry.instance().mapLayers()
for name, layer in layerMap.iteritems():
    if "point2line" in name:
        QgsMapLayerRegistry.instance().removeMapLayer(layer)

...but this only removed items from the layer panel not from memory.


texts: Adding a MapServer WFS layer in OpenLayers 3

texts:
I am following this example from OpenLayers' examples page so I can load a WFS layer pulled from MapServer. I copy and paste the code and the only thing I change is the EPSG in the function as the following snippet:

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(my_extent) {
      return 'http://localhost/cgi-bin/mapserv?map=/home/user/mapdata/mapFile.map&service=WFS&' +
          'version=1.1.0&request=GetFeature&typename=Rivers&' +
          'outputFormat=application/json&srsname=EPSG:2100&' +
          'bbox=' + my_extent.join(',') + ',EPSG:2100';
    },
    strategy: ol.loadingstrategy.bbox
  });

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

When I load it on the browser I get a white page and by opening the Network tab I click on my WFS link and a new page opens giving me the following error:

msWFSGetFeature(): WFS server error. 'application/json' is not a permitted output format for layer 'Rivers', review wfs_getfeature_formatlist setting.

I presume it is most likely a MapServer associated error but can somebody point me towards a solution?


texts: Google earth engine: filter a FeatureCollection by system id

texts:
I have a FeatureCollection in Google Earth Engine, where each Feature is referenced by an id tag. How do I filter the FeatureCollection based on that id?

Following attempts to filter id=3 did not work:

var g200 = ee.FeatureCollection('ft:11SfWB6oBS1iWGiQxEOqF_wUgBJL7Bux-pWU-mqd5')

print(g200.limit(5))

print(g200.filter(ee.Filter.eq("id", 3)))
print(g200.filter(ee.Filter.eq("system:id", 3)))


texts: How to run gdal_rasterize with sql option?

texts:
I'm running gdal_rasterize via a python pipe. When using the -l option, the raster is created just fine but I can't manage to get -sql option to work. Well my raster will be created but consists only of NAN's. What am I missing?

    rasterize_args = ['gdal_rasterize',
                      '-a', attribute_name,
                      '-sql', "select {attribute} FROM {layer}".format(attribute=attribute_name, layer=layer_name),
                      #'-l', layer_name,
                      '-a_nodata', str(no_data_value),
                      '-co', "COMPRESS=DEFLATE",
                      '-a_srs', srs.ExportToWkt(),
                      '-te', str(minx), str(miny), str(maxx), str(maxy),
                      '-ts', str(ds.RasterXSize), str(ds.RasterYSize),
                      ]
    if where_clause:
        rasterize_args.extend(['-where', where_clause])

    rasterize_args.extend([input_source, destination_raster])
    logger.debug("The gdal_rasterize cmd args are: %s", rasterize_args)
    ps = subprocess.Popen(
        args=rasterize_args,
        stdout=subprocess.PIPE
    )
    output = ps.communicate()[0]


texts: Is Labeling toolbar no longer available in QGIS 3.8.2 running on Mint 19.2 Tina?

texts: It seems the Labeling toolbar is no longer available with this version on Mint 19.2 Tina, though I have seen it on the Windows version of QGIS 3.8. Moving labels is a vital part of a large research and cartography project needed soon by my publisher.


texts: Why does field calculator compute area in square meters?

texts: I am using EPSG:102748, and have set my Settings/Option preferred measurement units to ft.  My CRS is NAD_1983_StatePlane_Washington_North_FIPS_4601_Feet, and when I use the measure tool, it returns the area in square ft or square miles or acres.  When I use the field calculator to compute an area field, I am just entering the $area function.  It is returning square meters instead of square feet.


texts: Reduce processing time for large GeoTIFF image

texts:
I am trying to process large GeoTIFF images (128 MB) using GeoTIFF reader provided by Geotools library. However the processing time is very high, can anyone suggest how to reduce this time.
Normal processing time for this image is 210000 ms, whereas my target is 500 ms.

GdalInfo of the image:

Driver: GTiff/GeoTIFF
Files: /Users/RGBI-7cm.tif
Size is 4000, 4000.
Coordinate System is `'
Metadata:
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4000.0)
Upper Right ( 4000.0,    0.0)
Lower Right ( 4000.0, 4000.0)
Center      ( 2000.0, 2000.0)
Band 1 Block=4000x1 Type=UInt16, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=4000x1 Type=UInt16, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=4000x1 Type=UInt16, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=4000x1 Type=UInt16, ColorInterp=Alpha

Code for processing the image in Java:

S3Object s3Object = awsFetchService.fetchObject(cropImageRequest.getAwsS3Bucket(),cropImageRequest.getPrimaryFilePath());
GeoTiffReader reader = new GeoTiffReader(s3Object.getObjectContent());
GridCoverage2D gridCoverage = reader.read(null);


texts: Populate a field from another field with a conditional

texts:
I am working in SDE geodatabase where I have a feature class 'Mange_Waste' and "Q1", "Q2" are two fields. I wanted to use python script for populating "Q2" from "Q1" with a condition. Like, if the value in "Q1" is "a", then it should return "1" in "Q2". I have hundreds of records in "Q1".

I am using following code however, it's giving me syntax error.

If Q1 == "a":
   Q2 = 1
else :
   Q2 = 0

Kindly tell me what I am doing wrong and what should be the correct code.


texts: How can I add a grid (from projection MTM8) on my layers

texts:
I need to add a grid on a map I created with my own layers and GPS data.

How can i do this?

I know how to do it in Arcgis, but not in QGIS.


texts: MySQL ST_Contains returning empty recordset?

texts:
Using MySQL 5.6, I have converted lat and lng fields into Points - using AsText, I can see this is OK.

I have another table to hold the polygon data, and again using AsText, I can see the correct points are in database.

My Query is this

SELECT P.PointId, P.Title, PLY.Name FROM polys as PLY,
points AS P
WHERE ST_CONTAINS(PLY.Polygon, P.gPoint)

There are points in the database within this polygon, no error, but none are coming through when I run this query.

Is this query correct, or do I need to use extra MySQL Extensions?


texts: Interpreting total cost variable in QNEAT 3 OD Matrix using QGIS

texts:
I am using QGIS 3.16.
The Matrix itself is optimised for Travel Time using pre-defined speed-limits on roads (These are ok, I have checked them). However, the total cost output is in the thousands (i.e., 2896) for some journeys, which makes no sense (if the results are to be interpreted in minutes) because the maximum distance away from the destination points is already 1 hour (a layer generated via ORS plug-in).
How should I interpret these variables, if not in minutes?


texts: How to open a Spatialite non-spatial view in QGIS

texts:
I'm using Qgis 1.8.0, and spatialite 2.3 in a Windows Vista machine.

I have a database in Spatialite, and I have created a view from one of its tables. When I open it in QGIS, using add Spatialite Layer (and enabling non-spatial tables), the attribute table shows ERROR in all of its lines and rows, but the number of rows and attributes seems correct.

I have tried both with spatial and non-spatial views withou sucess and opening tables (spatial or not) works well.

My table and view is something like this:

CREATE TABLE "Accoes"(
pkuid integer primary key autoincrement,
"UnidadeID" text,
"Tipo" text,
"Descricao" text,
"Data" text,
"Responsavel" text)

CREATE VIEW AccoesPoltest AS
SELECT c.* FROM Accoes as c

I have already read all info provided in this related question:
Can Qgis Read Spatialite Views?

I also found a very old bug ticket from QGIS, with the exact same issue, but it was closed as invalid: http://trac.osgeo.org/qgis/ticket/2178

This is one of my first approaches with Qgis + Spatialite, so probably something is missing.
Any idea? Thanks

Alexandre Neto


texts: Finding forest regions with consistently low value of NDVI

texts:
Using the code below, I have filtered my image collection(Landsat images) over the forest regions of India. (Forest regions were found using Hansen image). I have one mean image per month over 6 years from 2013 - 2018. I need to find out and demarcate those regions which have consistently low NDVI value (NDVI<0.25).

Since there are seasonal variations in NDVI value, a pixel corresponding to a particular point can have NDVI value greater or less than 0.25 at different times. So when I give a condition like updateMask(image.lte(0.25)) over a collection of images and then visualizing the regions on the map, I get all those regions where ndvi was less than 0.25 at least once. But I need only those regions where ndvi is less than 0.25 all the time.

What change should I make in my code?

// FILTERING OUT FOREST COVER IN INDIA

var treeCanopyCover = hansen.clip(india).select('treecover2000');
var greater25 = treeCanopyCover.gte(25);

// NDVI FUNCTION

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

// FILTERING THE COLLECTION

var startDate = ee.Date('2013-01-01');
var endDate = ee.Date('2018-12-31');
var filtered = l8.filterBounds(india)
      .filterDate(startDate, endDate)
      .filterMetadata('CLOUD_COVER','less_than',25)
      .map(addNDVI);

var maskedCollection = filtered.map(function(image){
  var x = image.updateMask(greater25);
  return x;
});

// FILTERING MONTHLY WISE OVER YEARS

var map_m = function(y)
{
  y = ee.Number(y);
  var months = ee.List.sequence(1, 12);
  var filtered_col = months.map(function(m)
{
var filt = maskedCollection.filter(ee.Filter.calendarRange(y, y, 'year'))
                .filter(ee.Filter.calendarRange(m, m, 'month'))
                .mean();

return filt.set('year', y)
      .set('month', m)
      .set('system:time_start',ee.Date.fromYMD(y,m,1)) ;
  });
  return filtered_col;
};

var years = ee.List.sequence(2013, 2018);
var monthlyimages = years.map(map_m).flatten();
var z = ee.ImageCollection(monthlyimages);

// // DISCARDING NULL IMAGES //

var nullimages = z.map(function(image)
{
  return image.set('count', image.bandNames().length());
})
.filter(ee.Filter.eq('count', 13));

// Finding regions where ndvi is less than 0.25

var ndviViz = {min: 0, max: 0.25, palette: ['red',  'blue']};
var ndviMasked = nullimages.map(function(image){
  image = image.select("NDVI");
  var x = image.updateMask(image.lte(0.25));
  return x;
});
Map.addLayer(ndviMasked.mean(), ndviViz);

I also have a doubt that whether it is right to use ndviMasked.mean() or ndviMasked.mosaic() while adding layer to the map.


texts: Getting Envelope Values From Extent using ArcPy?

texts:
Using following ArcPy snippet, I am able to get Extent of current map

import arcpy
import arcpy.mapping

## Sets the MXD file
IMXD = arcpy.mapping.MapDocument("CURRENT")
## Sets the Dataframe
DF = arcpy.mapping.ListDataFrames(IMXD, "Layers")[0]

print DF.extent

as:

  417939.472327849 241487.977474999 477757.540444938 363798.333225001 NaN NaN NaN NaN

can you please let me know how I can grab the Envelop value from This?


texts:
I'm trying to cut a line at the intersection, if a line crosses another line within a single multilinestring then and the end of the line is below a maximum distance threshold then it swaps the end of the line co ordinates to the intersection. As shown in this image.

The code I have so far is:

def trim_lines_to_intersection(lines, MAX_DIST):
'''If a road intersects with another road, end the road at the intersection
if the distance between the endpoint of the line and the intersection is less than
MAX_DIST'''

# initialize trimmed lines with list of original lines
#trimmed_lines = [line for line in lines]
trimmed_lines = [LineString(array) for array in lines]
#trimming_points = vertices_from_lines(trimmed_lines)
for i, line in enumerate(trimmed_lines):
    for q in [-1]:
        endpoint = Point(line.coords[q])

for road1, road2 in itertools.combinations(trimmed_lines, 2):
    x = road1.intersection(road2)
    if x.distance(endpoint) < MAX_DIST:
        Point(endpoint.coords[q]) = x
return trimmed_lines

Testing:

coords_SE = [((5, 5), (20, 5)), ((10, 10), (10, 0)), ((15,0), (15,10))]
lines_SE = MultiLineString(coords_SE)
test  = trim_lines_to_intersection(lines_SE, 6)

The error I am receiving is: SyntaxError: can't assign to function call.


texts: Extract Raster from Raster using Polygon shapefile in R

texts:
I am new to R and using the raster package. I have a problem extracting polygons from an existing raster file. If I use

extract(raster, poly_shape)

function on the raster it always creates a list with the data. What I really want is to extract another raster file that I can load with ArcGIS again.
After reading a bit more I think the crop function is what I really need. But when I try to use this function

crop(raster, poly_shape)

I get this Error:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

The files raster and poly_shape are the same for both functions. Can you tell me what could be wrong here? Is it even right that the crop function creates another raster and not a list?

EDIT: The extent() function does not work for me. I still get the same error. But I am sure the 2 datasets overlap! With the

extract(raster, poly_shape)

I get the right data from it. Just as a list and not as a raster like I want to have it. I just loaded the datasets in ArcGIS before and they fit very well so I did not check the projection. Now I tried

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

and you can see that the projections do not fit. The extract function seems to be able to automatically transform the files in the right way. I know that because I did the following:

I cut out the exact part of the polygon I extracted in R also in ArcGIS
I calculated the sum of all values of the extracted R polygon (list)
I calculated the sum of all raster cells that I cut out in ArcGIS

The 2 have the exact same result so I guess the conclusion should be that the extract function did work correct.
Now I have 2 options I guess:

I need a way to get a Raster out of the extracted list again or
The 2 datasets (raster + poly_shape) need to use the same prjection and the crop function should work

What would you suggest to do here?


texts: Optimising use of st_within for parallelisation

texts:
I'm trying to use st_within() to calculate points within polygons on a fairly complex shapefile. This operation tends to take 4-5 hours which is a bummer for debugging the code. I have access to a supercomputing cluster (120+ cores) and it seems sensible to make use of this tool to address this issue. However, it's unclear to me how I would even go about optimising these geospatial operations (from sf()) for multicore computing.

To get a sense of what I mean, you can try analysis using the UK National Forest Inventory: https://opendata.arcgis.com/datasets/bcd6742a2add4b68962aec073ab44138_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D

if (file.exists("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") == FALSE) {
download.file("https://opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D",
              destfile = "data/National_Forest_Inventory_Woodland_Scotland_2017.zip")
unzip("data/National_Forest_Inventory_Woodland_Scotland_2017.zip", exdir = "data")
}
forestinv <- st_read("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") %>% st_transform(paste0("+init=epsg:",27700))

You can see a concrete example in code here (especially lines 507-510)

Here's the key question-is it possible to parallelise these kinds of geospatial operations in R such as st_within? Even with a simplified set of geometries this process takes ages, but there aren't obvious ways to loop the operation and split across cores.


texts: How to create cohesive Spatial Pixels from Spatial Points Dataset

texts:
I have a Spatial Point DF spo (covering an irregular shaped area of interest). The data are not on a regular grid due to CRS transformation.

My goal is a raster with predefined resolution and extent of the area of interest (more spatial point data are to be mapped on this master raster).

Problems start when I

rasterize(spo, raster(ncol, nrow, extent, crs), spo$param)

I need to adjust nrowand ncol in a way so that I won't get moire patterns of NAs within my area of interest. I can't use a predefined (higher) resolution.

As a solution to this, I thought I would need some kind of Spatial Pixel DF spi, that covers my whole area of interest (just like meuse.grid in library(raster); data(meuse.grid)), and serves as a master grid. Then, I can use it to interpolate my data, e.g.

idw(param~1,spo,spi)

and by this, get full cover of my area of interest at my chosen resolution. But how can a SpatialPixelsDataFrame be produced from the point data?

So in my view, the question boils down to: How to produce meuse.grid from meuse dataset?

Maybe I'm taking the wrong approach here, so please let me know if more easily can achieved what I'm after, using a different way.


texts: GeoServer - change native SRS of cascaded WMS (other from Capabilities)

texts:
Short story:
How to change native srs of cascaded WMS service in GeoServer admin GUI? How to change native SRS from epsg 32633 to 4326? This native srs is blocked to change (it is somehow set automatically).

Long story:
I cascade WMS service with my GeoServer. Problem is that layer inside GeoServer sets "Native reference system" for some reason as EPSG:32633 (it is almost last in Capabilities, why it doesn't take epsg:4326, which is first???).

From this epsg:32633 I need to transform WMS to epsg:4326, but reprojection from native to declared gives me wrong BBOX coordinates, which ends in shift of my layer compared to the base map.

32633 BBOX coordinates are reprojected to 4326 are different, than those in Capabilities... Changing to other native srs can maybe force GeoServer to take good coordinates as they are in Capabilities.

<westBoundLongitude> 11.349143 </westBoundLongitude>
<eastBoundLongitude> 19.062675 </eastBoundLongitude>
<southBoundLatitude> 48.198780 </southBoundLatitude>
<northBoundLatitude> 51.668105 </northBoundLatitude>
<BoundingBox CRS="CRS:84" minx="11.349143" miny="48.198780" maxx="19.062675" maxy="51.668105"/>
<BoundingBox CRS="EPSG:4326" minx="48.198780" miny="11.349143" maxx="51.668105" maxy="19.062675"/>
<BoundingBox CRS="EPSG:32633" minx="244254.160680" miny="5342402.535258" maxx="798540.723911" maxy="5730646.738059"/>

When I tried to load WMS to QGIS I got the same "layer shift" error except when I set WMS projection to 5514 (this one made layer show correctly), which is not supported by GeoServer.

Here is list of supported epsg's written in GetCapabilities:

<CRS>CRS:84</CRS> <CRS>EPSG:4326</CRS> <CRS>EPSG:102067</CRS>
<CRS>EPSG:5514</CRS> <CRS>EPSG:28403</CRS> <CRS>EPSG:4258</CRS>
<CRS>EPSG:3035</CRS> <CRS>EPSG:32633</CRS> <CRS>EPSG:102100</CRS>
<CRS>EPSG:102113</CRS>

Capabilities are here!


texts: How to create tiles using GDAL from georeferenced tiff in C++?

texts:
I used GDAL to georeference a static image. However, when I load it in the map canvas (QGIS API), the image gets blurry when I zoom in. I need to generate map tiles but I cannot figure out how I would do that in C++. I read about gdal2tiles but I am still confused about how I would use the resulting .pngs. Can someone point me to any kind of example code or material that would make things clearer?

#include <gdal_priv.h>
#include <cpl_conv.h>
#include <cpl_error.h>
#include <ogr_spatialref.h>

using namespace std;
#include <string>
#include <cmath>

#include <QString>
#include <QDebug>

typedef struct Bounds{
    Bounds() : north(0), east(0), south(0), west(0) {}
    double north;
    double east;
    double south;
    double west;
}Bounds;

void GeoRasterizer::geoReference(const QString &source,
                                 const Bounds &bd,
                                 int imgWd, int imgHt,
                                 const QString &target)
{
    OGRSpatialReference osr;
    CPLErr error;
    char *wktProj = (char*)malloc(sizeof(char) * 50);
    GDALAllRegister();

    GDALDataset *srcDS = (GDALDataset*)GDALOpen(source.toStdString().c_str(), GA_ReadOnly);
    CPLAssert( srcDS != NULL );

    qDebug() << "Source " << source << " opened" << endl;

    GDALDriver *dstDrv = GetGDALDriverManager()->GetDriverByName("GTiff");
    CPLAssert( dstDrv != NULL );

    qDebug() << "Driver opened!" << endl;

    GDALDataset *dstDS = dstDrv->CreateCopy(target.toStdString().c_str(), srcDS, 0, NULL, NULL, NULL);
    CPLAssert( dstDS != NULL );

    qDebug() << "Target " << target << " opened" << endl;

    //Set the map projection
    osr.importFromEPSG(4326);
    osr.exportToWkt(&wktProj);
    error = GDALSetProjection(dstDS, wktProj);
    CPLAssert( error != CE_None );

    qDebug() << "Projection set to " << wktProj << endl;

    //Set the map georeferencing
    double mapWidth = abs(bd.east - bd.west);
    double mapHeight = abs(bd.north - bd.south);
    double transform[6] = { bd.west, mapWidth / imgWd, 0, bd.north, 0, -mapHeight / imgHt };
    error = GDALSetGeoTransform(dstDS, transform);
    CPLAssert( error != CE_None );

    qDebug() << "Transform set!" << endl;

    GDALFlushCache(dstDS);
    GDALClose(dstDS);
    GDALClose(srcDS);
}


texts: OpenLayers and MapServer: GetFeatureInfo and JSON

texts:
I'm using OpenLayers 2.0 and MapServer.

I would query MapServer and I would that it return to me a JSON file:

On my page, I generated dynamic URL as:

MAP/MAP.map&REQUEST=GetFeatureInfo&SERVICE=WMS&VERSION=1.1.1&LAYERS=TANK&FORMAT=image/png&TRANSPARENT=TRUE&SRS=EPSG:4326&BBOX=15.256938,38.195583,15.293846,38.210523&WIDTH=1720&HEIGHT=886&QUERY_LAYERS=TANK&INFO_FORMAT=text/plain&X=762&Y=685

In the file MAP.map (my Mapserver) I defined:

METADATA section:

  "wms_feature_info_mime_type" "text/plain"
  "wms_feature_info_mime_type" "text/html"
  "wms_feature_info_mime_type" "application/json"
  "wms_feature_info_mime_type" "application/vnd.ogc.gml"
  "wms_feature_info_mime_type" "gml"

LAYER section:
       HEADER   "/template/header.ftl"

[

TEMPLATE "/template/content.ftl"

{
    id: [TANK_ID],
    attribute: [DESCRIPTION]
}

FOOTER   "/template/footer.ftl"

]

If I launch the URL in my browser:

GetFeatureInfo results:

Layer 'TANK'
  Feature 100:

If I convert INFO_FORMAT in text/html, it works fine.

Why?


texts: Create animation of NDVI from ImageCollection in GEE

texts:
I'm trying to make an animation of NDVI images from sentinel collection in order to see clearly variance in the field.
The first part of the code works with little error: "Layer name is not a string: null",  but I want to use the second part of the code that is much better to work with.

/*  Copyright (c) 2018 Gennadii Donchyts. All rights reserved.
This work is licensed under the terms of the MIT license.
For a copy, see <https://opensource.org/licenses/MIT>.  */

var img1 = ee.Image('COPERNICUS/S2/20180824T081601_20180824T082948_T36SXA');
var img2 = ee.Image('COPERNICUS/S2/20180819T081559_20180819T082121_T36SXA');
var img3 = ee.Image('COPERNICUS/S2/20180814T081601_20180814T082955_T36SXA');
var img1_ndvi = img1.normalizedDifference(['B8','B4']);
var img2_ndvi = img2.normalizedDifference(['B8','B4']);
var img3_ndvi = img3.normalizedDifference(['B8','B4']);
var IMAGES = [img1_ndvi, img2_ndvi, img3_ndvi];
var animation = require('users/gena/packages:animation')
animation.animate(IMAGES, {label: 'label', maxFrames: 10})

Up to here the code works but that's a long way to do it, I would prefer to do it this way and add also visualization parameters

var visNDVI = {min: 0.3, max: 0.5, palette: 'FFFFFF, CE7E45, DF923D, F1B555'     +
'FCD163, 99B718, 74A901, 66A000, 529400, 3E8601, 207401, 056201, 004C00' +
'023B01, 012E01, 011D01, 011301'};

var geometry = ee.Geometry.Polygon(
  [[[-93.33501647520507, 38.46011907110406],[-93.33544562864745,     38.45339786087],
  [-93.33145450163329, 38.4532970379497],[-93.33117555189574,     38.45995104847996]]]);
var aoi =  geometry;
Map.addLayer(aoi, 'geodesic polygon');

var images = ee.ImageCollection("COPERNICUS/S2")
    .filterBounds(aoi)
    .filterDate('2019-01-01','2019-12-31');

var ndvi = images.map(function(image) {
  var result = image.normalizedDifference(['B8','B4']).rename("ndvi");
  return image.addBands(result);
});
imagesNDVI = ee.ImageCollection(ndvi);

var animation = require('users/gena/packages:animation')
animation.animate(imagesNDVI, {label: 'label', maxFrames: 10})


texts: Python error ("RecursionError") after updating QGIS

texts:
After updating from QGIS 3.12 Bucuresti to QGIS 3.14.16 with OSGeo4W, when opening QGIS (with or without a project), I get a Python warning "RecursionError" (sic!, QGIS says warning, Python says error). What does it mean and what do I need to do about it? I searched for solutions on the net, but I am less-than-rookie in Python, it gets technical  real quick, and I no longer see the forest for the trees.
WARNING    Traceback (most recent call last):
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 777, in _import
              new_path = os.environ['PATH']
              File "C:\OSGEO4~1\apps\Python37\lib\os.py", line 675, in __getitem__
              value = self._data[self.encodekey(key)]
              File "C:\OSGEO4~1\apps\Python37\lib\os.py", line 744, in encodekey
              return encode(key).upper()
             RecursionError: maximum recursion depth exceeded

Yes, it's only a warning, but warnings tend to come back later to bite us, so I'd like to not let it pass unnoticed.


texts: UnicodeEncodeError when creating Grass location on Windows

texts:
i´d tried to create a new location, but i have this in my console, iím learning to use grass gis, so can you help me?

   SIG\Archivos Grass\newLocation> no disponible

Iniciando GRASS GIS...
ERROR: No se encuentra la LocalizaciÃ³n
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS2~1.18\apps\grass\grass-7.2.2\gui\wxpython\gis_set.py", line 531, in OnWizard
    grassdatabase=self.tgisdbase.GetValue())
  File "C:\PROGRA~1\QGIS2~1.18\apps\grass\grass-7.2.2\gui\wxpython\location_wizard\wizard.py", line 2353, in init
    msg = self.OnWizFinished()
  File "C:\PROGRA~1\QGIS2~1.18\apps\grass\grass-7.2.2\gui\wxpython\location_wizard\wizard.py", line 2514, in OnWizFinished
    current_gdb = grass.gisenv()['GISDBASE'].decode(sys.stdin.encoding)
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\encodings\cp850.py", line 15, in decode
    return codecs.charmap_decode(input,errors,decoding_table)
UnicodeEncodeError: 'ascii' codec can't encode character u'\xed' in position 41: ordinal not in range(128)


texts: Variable road shield size determined from character length

texts:
In Maputnik / Mabox GL I am trying to display a shield around a road number, with the size of the shield determined from the length of the road number.
I am basing this on a working example as follows which looks up an image name e.g. "default_3":
["literal", "default_{ref_length}"]
This expression combines the word "default_" and the value stored in a field for the length.
In my data, I haven't got a field storing the number of characters. I would therefore like to automatically determine the length of the road number and use this in the same way.
The closest that I have achieved in terms of not receiving json errors is the following:
["string","default_",["length","ROAD_NO"]]
However, this doesn't produce the desired result. I suspect that this is because a literal is expected. I'm not sure how to present the expression as a literal though.
I know that a number is returned, because the following appears to return the correct number (I think!):
["length","ROAD_NO"]


texts: Javascript library for converting MGRS to Lat/Lng

texts:
Is there a JS library that converts MGRS to LL? I have found one that converts Lat/Lng to MGRS, but nothing that converts MGRS to Lat/Lng.

Someone asked this question in the past, but the answer was USNG to LL but I HAVE to do MGRS to LL.


texts: How to calculate area percentage from polygons from one shapefile within another?

texts:
I'm working on a project where I have to create python code that takes land-use data and adds a percentage field for each land-use category in a Census Tract shapefile. I've been trying to work with the steps manually in ArcGIS before developing the code, but I've had problems getting there.

Basically, I have a land-use shapefile with over 100,000 parcels categorized by land-use category (i.e., Residential, Commercial, etc.) as well as a Census tract shapefile. Within the Census Tract shapefile, I need to add a field with the percentage of area of each land-use category to be applied to each individual Census tract.

I've been working on this for hours and haven't gotten far at all. Could anyone offer any insight as to some steps I should try? I've tried intersecting, dissolving, joining tables (does not work) and more to no avail.


texts: trying to find areas within a certain driving distance along a road network WITHOUT a specific destination

texts:
I'm trying to use Arcmap 10.1 to find all the areas along a road network (a county roads layer in this case) that are a certain driving distance from a single origin point. In essence, this would be akin to a service area analysis, except there is only a single facility used rather than a whole list.

Applicability: trying to find where a subject could move to within the county (using census block layers) if they want to live no more than, for instance, an hour's drive from work/school/etc.

Any suggestions?


texts: Vector Legend item names in QGIS Map Composer

texts:
Is it possible to get the text for a legend item to go over multiple lines? I have some long item names and I have a set legend width that I have to use.

Without being able to make the text go over multiple lines I am unable to get my legend to fit the page properly.

Edit:

I have added an image to better illustrate the problem. As you can see the text for area type one goes out of the box I have set up for a legend. Personally I think short names in legends work best, but project managers do not always agree with this.


texts: Changing format of dates in ArcGIS 10.2 table

texts:
I am working with a table in ArcGIS 10.2 with over 500 entries containing dates. I am looking for a fast solution to transform the format of dates, from dd/mm/yyyy to mm/dd/yyyy.

So if a date is 10/08/1999, I need it to be 08/10/1999.

Is there any way to do this, maybe using the field calculator?


texts: Different number of available images between Scenes Explorer, Planet Explorer, and API

texts:
I'm attempting to download all PlanetScope 3 Band imagery that overlays ontop of an area of interest defined by a geojson file. The geojson file I'm using is seen here: https://gist.github.com/anonymous/621a279e1d6a56aa19844e9c99eb1e1d

When I go to the Scene Explorer and upload my geojson file as the area of interest, then adjust the timeline to "All Time", it says there are 36 scenes avialable.

When I go to the Planet Explorer and upload my geojson file as the area of interest, then adjust the timeline to "All Time", and set the type to types to "3 Band PlanetScope scene", it says there are 47 scenes available.

I have also used the API to conduct the same search. I created an API search filter using the geometry from the geojson file: https://gist.github.com/anonymous/f1943d3a1e9dbe5d4f7c2aaf631bcb1f . I conducted the search by using the 'download.py' script availbe through the Planet Docs. The code I used to conduct the search was:

download.py --query geom_filter.json --search PSScene3Band visual

Running this displayed that there was 64 available images.

So my questions are, why am I seeing different numbers for each search method?

And which method is delivering the accurate number of images?

Am I using any of the search methods incorrectly?

(Note: Searches were done on 2/1/2017)


texts: How to save the QGIS SCRIPT result in a folder

texts:
I would like to save my result from QGIS Script in a folder but I don't know how.

##Input_raster=raster
##result_out=output raster

outputs_SAGARASTERCALCULATOR_1=processing.runalg('saga:rastercalculator', Input_raster,[],'ifelse(gt(a,-1),1,0)',False,7,result_out)

I would like to save the results in a folder that I can easily change.

I am using QGIS 2.14.


texts: PostGIS Error: type "geography" does not exist

texts:
I just upgraded my development machine to Windows 7 and am therefore reinstalling PostgreSQL and PostGIS.  I used the Postgres Plus installer, as recommended on the Postgres website, and followed this tutorial for installation.  I created a database using the template_postgis and the application that I am currently developing was able to connect to the database with no problems.

However, my old code is not able to create tables using the PostGIS Geography like it used to.  For example, when I try to create the following table:

CREATE TABLE test_area (
    id SERIAL PRIMARY KEY,
    name VARCHAR(1000),
    area GEOGRAPHY(POLYGON, 4326),
    start_time DATE,
    end_time DATE
)

I get the following error:

ERROR:  type "geography" does not exist
LINE 4:         area GEOGRAPHY(POLYGON, 4326),
                     ^

If I look under the database in the Object brower of pgAdmin III, I see a bunch of Geometry functions, so the installation seems to have run correctly.  I'm not sure what my problem might be.  Is it possible for the installer to only install PostGIS Geometry, but not Geography?


texts: Is there a way of labeling the number of points with the same coordinates?

texts:
After reading through the following question (which is kinda related but not quite what I need):

How can I show points on the same coordinate as e.g. a bigger symbol in QGIS?

I found out that there might be a need to label these locations that have multiple points with same lat/long values. For example, if a location has 20 features in it, I want the label to say 20. Or even more useful: if a location has 20 features in the location, and 15 fulfill a condition (string = abc), I want the label to say 15.

Is there a way of doing so in QGIS?


texts: Why do I get the same error code when Geocoding?

texts:
Whenever I use the 'MMQGIS' plug in to 'Geocode' a CSV file containing postcodes to get the longitude and latitude which will enable me to plot the locations as a layer on QGIS, it always come up with the error code seen below.

Failure creating output file: Creation of data source failed (OGR error: C:\Users\louis.gorringe\Desktop\Chemical.shp is not a directory.)

So although the shapefile may be present after running the geocoding. whenever I exit QGIS versions 3.10.0 and re open it, all shapefiles have dissapered. Why is this happening?


texts: Unable to perform SAGA point snapping to line on Win and Mac

texts:
I have been trying to snap points from a shapefile to a line feature from another shapefile. However, each time I do this, I get the following error(s) (depending upon the OS used):

(For OSX El Capitan):

File"/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/saga/SagaAlgorithm.py", line 341, in processAlgorithm
              with open(prjFile, 'w') as f:
             PermissionError: [Errno 13] Permission denied: '.prj'

The above error being common to Mac and Windows as well. Sometimes while running the tool on windows environment is able to pass the permission error stage (somehow and its pretty weird, only sometimes) and generates only '.prj' file but fails stating:

  'not able to generate shapefile'

I have tried it on QGIS Versions: 2.18, 3.2, 3.4.3 and 3.4.4.

All the supporting packages (GDAL) and respective python versions have been properly installed (3.6.7 for QGIS 3.0+).

My system specs are as follows:

System: Macbook Pro Mid 2012 i5 El Capitan + Windows 10 Pro (Bootcamp)

How can I fix this error?


texts:
I want to exclude cells with 3 specific values, what's the proper syntax?

I've got "Attribute" != 'Value1', which works fine.  But what's the right syntax to add the other values?  I thought I could just use OR followed by the next value, but that doesn't work, e.g. "Attribute" != 'Value1' OR 'Value2' but that isn't working.


texts: QGIS 2.18 can't read multiple rasters in on geopackage

texts:
I'm trying to add multiple rasters in one geopackage

But QGis 2.18 seems to be unable to deal with...

To add the rasters I use GDAL commands :

gdal_translate -projwin 578000 6270000 580000 6268000 -a_srs EPSG:2154 -of GPKG -co RASTER_TABLE=image1 -co TILE_FORMAT=JPEG d:\donnees\ortho.ecw c:\temp\test.gpkg
gdal_translate -projwin 579000 6275000 581000 6273000 -a_srs EPSG:2154 -of GPKG -co RASTER_TABLE=image2 -co APPEND_SUBDATASET=YES -co TILE_FORMAT=JPEG d:\donnees\ortho.ecw c:\temp\test.gpkg

Everything seems ok

But when I try to open the raster from QGis, I get this dialog but then nothing happens...

If I add only each image in one geopackage, it is ok but I would like to put 2 images in one geopackage.


texts: Prevent labels being cut off / clipped at edge of map

texts:
I've read a few other issues about trying to prevent labels from rendering if they run off the edge of the map canvas. This is particularly a problem in Composer, where I get lots of country labels which are only half on the page.

https://issues.qgis.org/issues/10314

In the above post in the QGIS issue queue, Nyall Dawson says "There's an option in the global label settings to prevent clipped labels at the edge of maps".

I can't find this setting. Does anyone know where it is? Or any other solution for this problem?


texts: MRT (MODIS Reprojection Tool) input error

texts:
This is my first time to use Modis Swath Tool to deal with images. I've downloaded 4 MOD44B granules (using the FTP script method), and am trying to use MRT to mosaic them. When I add the 2 files I get an error saying firstly One of the input files does not exist or the inputs are not the same data product type then Input file does not exist!. If I add files individually, I get an error saying Error in Module: ReadHrFile, Opening Input Header File.

I saw this error before in this website, but I still failed to fix it following the answer to this question by opening ModisTool.BAT rather than the ModisTool.JAR executable. Any one knows how to fix it? Any help will be appreciated. Thanks very much.

Versions:

MRT: 4.1 (March 2011)
Windows: 7
Java: 7


texts: Editing feature geometry in action using PyQGIS

texts:
I have fieldwork photos stored in a vector layer, geotagged from EXIF GPS data. Sometimes the GPS is off, or even missing, and I need to move their position manually. Unfortunately, the Advanced Digitization Move feature command requires I click on the current (geotagged) position (after selecting the feature first) and only then the corrected position, which is awkward if the geotagged position is say Point (0 0).

As a workaround, I use the Lat Lon tools plugin to copy the desired position to the clipboard (as WKT), and then I've written a PyQGIS feature action to paste that as the feature's geometry. However, following the PyQGIS cookbook the code below uses the dataProvider() to change the geometry, and so the currently open form and attribute table that are invoking the action don't pick up on the change. Even hitting the refresh icon isn't enough.

How could the action code be changed to use the (already active) QGIS editing session instead of going behind its back to the provider?

# find ourselves, layer first
ourLayerId = '[% @layer_id %]'
ourLayer = QgsProject.instance().mapLayer(ourLayerId)

# now the feature, note getFeatures returns an iterator, but will be unique
fid = [% $id %]
featIterator = ourLayer.getFeatures(QgsFeatureRequest().setFilterFid(fid))
feature=featIterator.__next__()

# figure out new geometry, using WKT from clipboard
# should check CRS (if WKT is output of Lat Long tools it should be 4326!). For now assume it's fine.
geomWkt = QgsApplication.clipboard().text()
geom = QgsGeometry.fromWkt(geomWkt)

# make the change. This bypasses the editing cache going to the data provider
result = ourLayer.dataProvider().changeGeometryValues({fid : geom})


texts: Calculating distance to nearest object using PostGIS?

texts:
I would like to calculate the shortest distance between each building and nearest sewage pipe. I followed steps described on CARTO blog but I encounerd an error.
I have data in same database but in different schemas:

table name for buildings is ks2_st and schema is public
table name for sewage pipes is t2100_ill and schema is gji

The query is:

SELECT
  a.mid_st,
  a.geometry,
  ST_Distance(geography(b.geometry), geography(a.geometry)) as distance
FROM
SELECT DISTINCT ON (a.geometry)
  FROM public.ks2_st a) as buildings
CROSS JOIN LATERAL
(SELECT b._mid_gji, b.geometry
  FROM gji.t2100_ill b
  ORDER BY a.geometry <-> b.geometry LIMIT 1) as sewery

When I try to run query I get message:

  ERROR:  syntax error at or near "FROM" LINE 7: FROM public.ks2_st a)
  as buildings

Can you please help me find the error, because I do not know what the error could be?


texts: GTOPO30 elevation value are wrong in QGIS

texts:
I downloaded GTOPO30 DEM archives from USGS(earthexplorer.usgs.gov/), and I just simply drag the file into QGIS, but somehow I find the elevation isn't the same with the actual scenery.
I want to draw the terrain profile but the value isn't correct.
I open the properties→ Metadata in QGIS and the maximum value is different from the official document.

Example: GTOPO30 w180 n90

The official specification is attached.
The maximum value of this dem should be 6098, but at Qgis the Max is 4659.

Could anyone help me to deal with this error?


texts: CartoDB import api: upload data from stdin instead of file

texts:
I am trying to upload data using the CartoDB import api. Uploading a regular file works but my attempts to upload data from stdin are failing.

The example of uploading a file works perfect.

curl -F file=@/Users/mbj/Downloads/cvr_count.csv "https://{CARTODBUSER}.cartodb.com/api/v1/imports/?api_key={APIKEY}"

I want to load data comming from stdin instead.

According to curl man page

  To read content from stdin instead of a file, use - as the filename. This goes for >both @ and < constructs.

To test this:

cat /Users/mbj/Downloads/cvr_count.csv | curl -F file=@- "https://{CARTODBUSER}.cartodb.com/api/v1/imports/?api_key={APIKEY}"

Response:

{"item_queue_id":"7b8477f9-0926-4f11-b898-c9d9f149ae00","success":true}

But the table never appears in CartoDB. Has anybody tried this?
Am I doing something wrong ?


texts: Can I collect Seismic data with an array of cell phones?

texts:
Suppose I have an array of cell phones, each having:

an accelerometer
a clock
a GPS receiver
Windows Phone 7 (WP7)

How hard would it be to write a program so that synchronizes the phones to a common clock (preferably a GPS clock)?

How hard would it be to write a program that starts recording accelerometer readings when it receives text message?  I would want it to collect for a duration specified in the text message, and then forward the collected data (including location) back to an address specified in the text message.

How hard would it be to write a program (on a PC or a cell phone) that recieves the data, processes it to compute a "fix" (ellipse) on the source signal location, then sends that fix back to each cell phone as a text message?

Alternatively, could I record audio spectrograms instead of accelerometer readings?


texts: Changing reference scale of annotation feature class in ArcGIS Pro?

texts:
I've looked all around, but can't find a way to change the reference scale of the Annotation FC in ArcGIS Pro.

Is there a way to do it instead of having to create many FCs?


texts: How to create a ClassBreakRenderer without data coming from Feature Layer?

texts:
I would like to create a ClassBreakRenderer to represent integer values. The symbols I am using are defined as PictureMarkerSymbols with JSON and they are just a filled blue circle with different sizes that were created with ArcMap 10.1 and published into ArcGIS Server. This is just to introduce the thing.

The problem is that in all the examples I see on the ArcGIS Javascript API about Rendering, the examples do always get the data from a particular column of a Feature Layer, but my data do NOT come from this data source: they come from an array of integers in client side. Just the symbol comes from the server-side, but not the integer value.

Is there any way of "tricking" the constructor of ClassBreakRenderer in order to accept an array (or collection, dictionary, set, whatever) of integers to be classified in symbols? There should be "life" a part of Feature Layers, isn't it?

Thanks for your help!

EDIT

The thing is that my data source is an array of integers like this:

var myArray = [15, 37, 47, 12, 125, 87, 36, 66, 10, 4, 18, 19, 23]

AND NOT:

var featureLayer = new esri.layers.FeatureLayer("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/3", {
          mode: esri.layers.FeatureLayer.MODE_SNAPSHOT,
          outFields: ["*"],

And in principle, the ClassBreakRenderer receives a column of that Feature Layer as:

var renderer = new esri.renderer.ClassBreaksRenderer(symbol, "POP07_SQMI");

But I do NOT have that column "POP07_SQMI" I have an array called "myArray" and my question is that if it is possible to use the array instead of the column name that I don't have. Because world should not be limited to the existence of feature layers to retrieve data from them.... :S


texts: Calculate coordinate for a point along given line

texts:
Using Openlayers or plan Javascript, how would i get the coordinate of a point that is 500 meters along this line?
Im trying to calculate the coordinate (lat long) x meters along a given line

Example coordinate for a line :

[59.8736731,10.803402600000027]
[59.9319572,10.688514700000042]


texts: Creating a GeoExt Action from OpenLayers Function

texts:
I am trying to create a GeoExt function out of this OpenLayers function, which displays the getFeatureInfo.

I was able to create this GeoExt.Action() but it does not give me any output. Any idea where am I going wrong? This is the function addToPopup().

This is my application, http://128.196.142.12/geo/test/test_new.html.


texts: How can I get the FAO gaul Adm name with special characters?

texts:
I am building an App and I am calling FAO gaul from my backend (Python) to get the value like this:
return  (ee.FeatureCollection('FAO/GAUL_SIMPLIFIED_500m/2015/level0')
                .filter(ee.Filter.eq('ADM0_NAME',variable)))

The issue comes when I try to get the value of countries with special characters such as Côte d'Ivoire and instead of that value my backend gets C�te d'Ivoire.
Anyone knows why I get this question mark and how to get the actual character?


texts: Conversion from Namibian Coordinates to Google Maps (WGS84)

texts:
I have coordinates in Namibia, but don't know how to convert them to use them in Google Maps.

The system is Lo22/19°.

Then Y:  +55 435,89  X: -498 157,28

Can anyone help with a conversion?


texts: PostGIS: Geometry column registering as GeometryCollection, not converting to MULTIPOLYGON

texts:
I'm trying to aggregate data from a tabular data source and then join it to a geometry file (sl_adm2) in PostGIS using a primary key. The geometry file was imported using shp2psql.

CREATE TABLE example_geom AS
SELECT sl_adm2_sums_alias.* , adm2_code as district_pcode
FROM gis_schema.sl_adm2
JOIN (
    SELECT st_collect(sl_adm2.geom) AS geometry, sum(sl_sos.HEALTH) AS health_sum,
        sum(sl_sos.PROTECTION) AS protection_sum, sum(sl_sos.WASH) AS wash_sum,
        sum(sl_sos.cccm) AS cccm_sum, sum(sl_sos.food) AS food_sum,
        sum(sl_sos.livelihood) AS livelihood_sum,
        sum(sl_sos.education) AS education_sum,
        sum(sl_sos.agriculture) AS agriculture_sum,
        sum(sl_sos.HUMANITARIAN_ACCESS) AS humaccess_sum,
        sum(sl_sos.logistics) AS logistics_sum,
        sum(sl_sos.SHELTER_AND_NFI) AS shelter_and_wifi_sum, sl_sos.district
    FROM sl_sos
    JOIN gis_schema.sl_adm2 ON sl_sos.DISTRICT_PCODE = sl_adm2.adm2_code
    GROUP BY sl_sos.district)
sl_adm2_sums_alias ON
sl_adm2.adm2_name = sl_adm2_sums_alias.district;

and then

ALTER TABLE example_geom
ADD PRIMARY KEY (district_pcode);

The aggregation and join work fine but the geometry column in the resulting table does not appear to work. Further investigation reveals that the geometry column is registering as a GeometryCollection. I tried to change this with the following code.

ALTER TABLE example_geom ALTER COLUMN geometry
SET DATA TYPE geometry(MultiPolygon) USING ST_Multi(geometry);

This gave me the below error

ERROR: Geometry type (GeometryCollection) does not match column type (MultiPolygon)
SQL state: 22023

I've tried to fix this or get around it with no success. Could anyone help?

I'm using Postgres 9.4.1 with PostGIS version 2.1.5 on a Mac OS X 10.10.2.


texts:
I want to extract contour from scanned map, but that map is scattered very much, so it is impossible to extract automatically.

So I make line shapefile and try to use tracing tool

but I cannot do that because my map is just image.

How can I use tracing tool for my map?


texts: Error in change of projection from Lat/Long to UTM

texts: Why does the the DEM when converted to UTM Zone 43N system give all values as 0 in QGIS? I followed the standard procedure of saving the parent layer and then assigning the CRS. When I open the same in ENVI, after editing the header and reassigning the projection of the parent layer, I can see that the units are in metres, but the same when opened in QGIS renders geographic lat/lon system.


texts:
I am currently working with high volume ArcGIS (9.3.1 & 10) production of PDF files, part of the companies requirement is that a database (Access 2003) needs filling out with each map, this needs an automated (sequential) map number, the mxd location and the user (3 digit reference). This also needs to be filled on the map template in the specified fields of the map.

I have seen some python script for exporting the layers of a map to a table, is there a way that I can automate my process so that I can fill the "username" field on the map and then have ArcGIS populate the "figure number" field with the next sequential number and then populate a table which I can link to the DB which my company uses?

This would make my life a lot easier and I can see how it is possible though my python is letting me down.


texts: How to use create a plugin that contains a set of processing scripts using Qgis Script feature?

texts:
I have some script that I would like to share easily.

I have read here that its possible to create a plugin using a set of scripts using "Scripts/Tools group in the Processing toolbox, double-click on the Create script collection plugin item"

Unfortunately I don't have this tools in scripts in ubuntu 16.04 with QGIS 2.14.3, or in windows.

Does anyone know how to solve this problem because if I go into QGIS source I have seen this code that is missing in my QGIS installation: CreateScriptCollectionPluginAction.py.

Now this tools is available in Qgis and works well when you apply the following change (@Joseph helped):
C:\Program Files\QGIS 2.16.1\apps\qgis\python\plugins\processing\script
And then change the line after the imports:

WIDGET, BASE = uic.loadUiType(
    os.path.join(pluginPath, 'ui', 'DlgConfig.ui'))
with

WIDGET, BASE = uic.loadUiType(
    os.path.join(pluginPath, 'ui', 'scriptselector.ui'))
Re/start QGIS and hopefully it will work:

It works well because when I go to plugin I can enable it. However, unfortunately, after that I am not able to see any new menu where to find the script, whatever in Plugins, processing option.

Did you know where to find the menu ?


texts: QGIS Export: Set DXF Layers By Attribute

texts:
I have a shapefile layer of contour lines that I'd like to export to DXF, but I want to have each feature be on a layer defined by an elevation attribute. I saw here that something like this should be possible, but can't find how to do it.

I tried both "Save As" and "Project -> DXF Export," and saw no obvious option to set the DXF layers. I also tried creating an attribute called "layer" (and also tried "Layer" and "layers") to see if either export method would automatically use this attribute for the DXF layer, but got nothing.

I can get this to work using ogr2ogr, but is there a more direct way to set the DXF layers on export in QGIS? (BTW I am using QGIS 2.18.3)


texts:
I updated QGIS and accidentally rewrote 2.18 as ltr with 3.4 Madeira.
My apps work with older version. This means there must be qgis-full and qgis-ltr with version 2.18 to make them work correctly.
How can I restore my previous installation or find packages for installation from local folder.
I use Windows 10 64 bit.


texts: Zooming specific WMS using OpenLayers 2?

texts:
I'm looking for a way to use multiple WMS on my application. If a user zooms in/out another baselayer is supposed to be selected. As for now, I've come across the following:

map.events.on({ "zoomend": function (e) {
    if (this.getZoom() > 10) {
      mapnik.setVisibility(false);
      layertest.setVisibility(true);
    }
    else {
      layertest.setVisibility(false);
      mapnik.setVisibility(true);
         }
                                        }
            });

The problem I still encounter is really bad performance. Well it works, but is there maybe a faster way of doing this?


texts:
I would like to import a csv file to my Carto account using R. I know the csv has been generated, and I'm able to save it in other locations.

Running the code to upload this file to Carto appears to be successful (see success message below), yet I'm unable to find/retrieve the csv in my carto account.

Here is the code I'm using to import my csv:

apiurl <- "https://"MY ACCOUNT".carto.com/api/v1/imports/?api_key= "MY API KEY"

inputFile <- "MyFile.csv"

writeLines(inputFile)

resp<-POST(apiurl, body=upload_file("MyFile.csv"), encode="multipart")
content(resp)

I'm pretty sure the code I'm using is correct, since it generates this message:

content(resp)
$item_queue_id
[1] "(a unique alphanumeric string) "

$success
[1] TRUE

... but where is my csv?


texts:
This is related to this previous post.

As was trying to implement the said answer and inserting it to my code, but the code doesn't seem to work.

#importing all the necessary modules
import os, glob, gdal
from subprocess import call
#call(["ls", "-l"]) #calling the "ls -l" command in linux

shp = "/home/brent/Documents/Hansen_new_data/Palawan_only.shp"

os.chdir("/home/brent/Documents/Hansen_new_data/New_data/10_110/NAMRIA_based")
extent = ("116.928337 7.506939 120.000087 10.000189")
for rname in glob.glob("*AMRIA_gain_10N_110E.tif"): #iterate through the files using glob.glob
    raster = os.path.join("/home/brent/Documents/Hansen_new_data/New_data/10_110/NAMRIA_based", rname) #join the shp file name to a new file path name
    raster2 = rname.replace("NAMRIA_gain_10N_110E","Palawan")
    outRaster =  os.path.join("/home/brent/Documents/Hansen_new_data/New_data/10_110/Province/gain/", raster2) #specifying the output raster new name
    warp = 'gdalwarp -te \'%s\' -cutline \'%s\' -crop_to_cutline -overwrite -co compress=LZW \'%s\' \'%s\'' % (extent, shp, raster, outRaster) #calling gdal warp, crop to cutline is specified here
    os.system(warp) #running gdalwarp, uncomment this line to run


texts: Adding a field to a shapefile but won't let me type manually

texts: When trying to add a new field to a shapefile in ArcCatalog, I have been told to type the field name manually into the vacant column in the shapefile Properties> Fields tab. However, the system won't let me type in the row under the field name column. I type and nothing shows up.


texts: Classify a slope raster in Image Analysis

texts:
I have a DEM in a Mosaic Dataset covering a large area (300 x 200 km). I want to create a classified slope raster showing slopes of say 0 - 10 degrees, 10 - 20 degrees, and then 20 - 90 degrees. Because the dataset is so large I want to avoid having to create a separate slope raster dataset, but instead use the Image Analysis so I can visually see where the slope classes are.

In the Image Analysis window I can create a Slope Function for my dem and display it using the Stretched color ramp in the Symbology tab of the Layer Properties. I want to use the Classified renderer, but Arc requires the data to have a histogram. When I say Yes, it chugs away computing the statistics and histogram. For really large datasets it can take a very long time. I'm wondering where the histogram and statistics are stored for a Image Analysis layer because it this Compute Histogram dialog keeps re-appearing.

I thought there could be a way to create this with a Function chain. There are functions called: Classify, Colormap, Remap that sound promising, but I haven't gotten any of them to work

The Classify function wants a .ECD or .ACT file as input. How do I create a ECD or ACT file or is there a better way to accomplish this?

Function Chains:


texts:
I have a polygon layer and a point layer.
1) I wrote a spatial query to count how many points fall in each polygon:

SELECT polR.codeR, count(pntAZ.geom) as tot_pntAZ
FROM boundaryR as polR LEFT JOIN waypointsAzim as pntAZ
ON st_contains(polR.geom,pntAZ.geom)
GROUP BY polR.codeR;

This SELECT does the job correctly.

2) I want to create a view with the centroids of boundaryR polygon layer, attaching also the result of the previuos SELECT statement. To do this I wrote the following:

CREATE OR REPLACE VIEW view_waypAZ AS
SELECT polR.codeR, polR.nameR, ST_centroid(polR.geom) as geom
FROM boundaryR as polR;

but the result is a non-spatial view: it's just a table, but I need a spatial view of points (the centroids). So I tryed to add the geometry column:

INSERT INTO geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, "type")
SELECT '', 'geodata', 'view_waypAZ', 'geom', ST_CoordDim(geom), ST_SRID(geom), GeometryType(geom) FROM geodata.view_waypAZ LIMIT 1;

The query runs succesfully but the added row in geometry_columns table is not complete: SRID is 0 (should be 4326) and TYPE is "geometry" (should be MULTIPOINT).

So, the question is: how can I make a spatial view of centroids from a polygon layer? And how to "attach" to it also the results of statement 1? I hope everything is clear enough.
Thanks in advance


texts: Importing LiDAR point cloud classes in Potree?

texts:
I can't find a way to import classes (roof, vegetation, ground, etc..) in Potree.

I generated the web page, except I get the black point cloud. When I, for example, open .las file in ERDAS Imagine, I get coloured points.

I tried to type these commands:

PotreeConverter /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter -o /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter/output --output-attributes CLASSIFICATION --projection 4326 --generate-page MYPAGE

PotreeConverter /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter -o /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter/output -a CLASSIFICATION --projection 4326 --generate-page MYPAGE

PotreeConverter /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter -o /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter/output --output-attributes RGB CLASSIFICATION --projection 4326 --generate-page MYPAGE

PotreeConverter /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter -o /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter/output -a RGB CLASSIFICATION --projection 4326 --generate-page MYPAGE

PotreeConverter /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter -o /home/user/LAStools/LASzip/PotreeConverter/PotreeConverter/output --projection 4326 --output-format LAS --generate-page MYPAGE

But all I get is black cloud point. I also tried to colour by height, but it's also black.


texts:
How would you focus your ArcMAP application window using VBA?

I have a script that runs and selects polygons and zooms to them, but I would like ArcMap to focus and come to the front at that point.

Suggestions?

Dim pDoc As IDocument
Dim pApp As IApplication
Dim pAppROT As IAppROT

Set pAppROT = New AppROT

If pAppROT.Count > 0 Then
Dim p As Integer
Dim x As Integer
p = 0
x = 0
Do While p < pAppROT.Count
If TypeOf pAppROT.Item(p) Is IMxApplication And (pAppROT.Item(p).Document.Title Like "Report.mxd") Then

Set pDoc = pAppROT.Item(p).Document
Set pApp = pDoc.Parent
x = 1
pAppROT.Item(p).Visible = True

End Sub

Set focus = pAppROT.Item(p)

Exit Do
Else

End If

p = p + 1

Loop

End If

If x = 0 Then

    Set pDoc = New MxDocument
    Set pApp = pDoc.Parent

    pApp.OpenDocument ("c:\Report.mxd")

    pApp.Visible = True

End If'


texts:
I uploaded my own symbols to use in the legend of my map.

However, after I refresh the map some of the symbols will disappear. And each time after that I refresh the map, more symbols disappear. Right now I just have the symbols saved in a folder under "My Content" tab.

Do I need to use another online source in order for them to save correctly?


texts:
I tried to follow what was suggested in How to use OS Terrain 50 in QGIS, but I got the following problem:

in spite of being generally able to generate .vrt files using gdalbuildvrt from unzipped .asc files, when I try to do it using /vsizip/ the resulting .vrt seems corrupted;
it only contains NaN values, as can be seen from the following QGIS screenshot:

GDAL version: GDAL 1.11.2, released 2015/02/10

First few lines for a single .asc file input:

ncols 200
nrows 200
xllcorner 440000
yllcorner 1200000
cellsize 50
    -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9 -0.9

First few lines of the corresponding output:

<VRTDataset rasterXSize="200" rasterYSize="200">
  <SRS>PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000],PARAMETER["False_Northing",-100000],PARAMETER["Central_Meridian",-2],PARAMETER["Scale_Factor",0.999601272],PARAMETER["Latitude_Of_Origin",49],UNIT["Meter",1]]</SRS>
  <GeoTransform>  4.4000000000000000e+05,  5.0000000000000000e+01,  0.0000000000000000e+00,  1.2100000000000000e+06,  0.0000000000000000e+00, -5.0000000000000000e+01</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="0">/vsizip/hp40_OST50GRID_20150709.zip/HP40.asc</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="200" RasterYSize="200" DataType="Float32" BlockXSize="200" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="200" ySize="200" />
      <DstRect xOff="0" yOff="0" xSize="200" ySize="200" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Does anybody know what might cause it?


texts: Creating layer showing difference between two shapefiles using ArcMap

texts:
I have 2 zoning shapefiles for different years. I am creating a swipe map that will show those two years separately to visualize change in zoning. In addition to this, I want to create a third map that shows the difference between them.

Which tool or which approach would accomplish this analysis?


texts: Leaflet unproject from EPSG:25831 doesn't work

texts:
In a leaflet map that uses EPSG 25831 to get WMS request, it works fine to call to "project" function to get UTM coordinates of one coordinates pair that I got from the map in WGS84, but inverse function to "unproject" a pair of coordinates in UTM in order to get this equivalent in WGS84 doesn't work.

The console message is "unproject  is not a function"

The example code:

var crs = new L.Proj.CRS('EPSG:25831',
'+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
{
    resolutions: [2000, 1100, 550, 275, 100, 50, 25, 10, 5, 2, 1, 0.5, 0.25]//,
    //origin: [0, 0]
}),
map = new L.Map('map', {
    crs: crs,
    continuousWorld: true
});

L.tileLayer.wms("http://server/service?", {
layers: 'topographic',
format: 'image/png',
transparent: true,
continuousWorld: true
}).addTo(map);

map.setView([41.3, 1], 6);

// UTM to 4326
var point = L.point(4651000,400000);
var newPoint = crs.unproject(point);
console.info(newPoint);

In the other hand, "project" function, works fine:

var testPoint = new L.latLng(42.12, 1);
var geo_point = crs.project(testPoint);
console.info(geo_point);

Versions used are, Leaflet: 0.7.7, proj4leaflet: 0.7.2, proj4: I don't know, it was at the Proj4Leaflet-0.7.2.zip from: github.com/kartena/Proj4Leaflet/releases

Any suggestion?


texts: Leaflet draw: prevent drawing through a custom overlay

texts:
I have a leaflet draw implementation (in Angular) with a custom sidebar that overlaps the Map. It should not be possible to draw through that sidebar.
For circlemarkers I could solve this problem easily by using

pointer-events: all;

and

event.stopPropagation();

for the clickevents. That way any clicks will not arrive on the map, so no circlemarker is drawn.

My problem: this solution does not work when drawing eg. rectangles. Even with the code above, any click on the sidebar starts drawing a new rectangle or finishes the drawing of a rectangle.

How can I prevent the Map from registering these events on top of the sidebar. Why does it work for circlemarkers an not for polygons?


texts: Using gdalwarp for warping Space View Perspective or Orthographic tiff?

texts:
i need to warp an image with a Space View Perspective or Orthographic project. What is the parameter for the gdalwarp -s_srs "+proj=???? +ellps=WGS84 +datum=WGS84" -t_srs EPSG:4258 in.tiff out.tiff ?
I use +proj=longlat for Latitude/longitude (gridDefinitionTemplateNumber = 0).


texts: Tessellation in Teradata is deprecated, but Geoserver expects it to improve spatial query?

texts:
In the project am involved with, I am responsible to:

store geospatial data in Teradata database
join the geometry with other non-geospatial attributes in the database
publish them with GeoServer with WMS

All the above requirements are passed actually. However the performance of WMS is noticeably slow, much slower than PostGIS-backed Geoserver.

So I dug into GeoServer log to see what happen behind the screen during WMS request, and i found several lines similar to the below:

  2018-01-23 11:47:22,402 INFO [data.teradata] - Tessellation info not
  available for GEOM, unable to perform spatially indexed query

It seems that the performance is slow because i did not setup Tessellation in Teradata database.

Well, I did not setup Tessellation because it is deprecated. Teradata documentation titled Teradata Database SQL Geospatial Types 14.10, April 2014, (PDF) page 4, says that since the version 14.10 Native indexing (NUSIs) for geospatial data replaces tessellation as the preferred method for improving geospatial performance. Tessellation is deprecated.

So there's something funny here:
Teradata document says that Tessellation is not the prefered method for improving performance, therefore is deprecated; but GeoServer expects Tessellation for spatial index.

Can anyone explain this? Or do i simply set Tessellation anyway?

The ecosystem used is:

Teradata Database 16.00.00.14 on quite big machine.
GeoServer 2.10.1 on JBOSS EAP 6 on RedHat Enterprise Linux 6 64-bit


texts: Displaying raster image (tiff, img, hdf) using gdal, Tkinter?

texts:
I am building modis data processing tool in python. I am building GUI using Tkinter and installed various geospatial library. i have tried to visualize image using matplot lib but it is showing image in new window. i want make tool like QGIS, or GRASS.

How do they build raster viewer?

installed python 2.7.13, GDAL/OGR, Numpy, (all basic lib).


texts:
We have an internal GIS library which has a very wide range of data from numerous sources.  We use this to supplement information provided by a client.  Our internal GIS files have coverage for the country, however for specific projects I only need a small subsection of the data.

At the start of a project I like to consolidate all of the information togethor. Is there a method to quickly copy an extract of the data to the project directory - to minimise the data set and speed up the work and also to make a copy for the records.


texts: Trying to run GDALbuildVRT in command line

texts:
I have a directory of OS data stored as tiff files (approximately 2000), which I want to create a virtual raster with. I understand that you can run this in the QGIS UI but because of the large number of files it is not possible (see: QGIS virtual raster - insufficient permissions)

I therefore want to run this in the command line of which I have very little experience. My workflow so far is this:

I have changed the directory to the "data" folder where all of the raster tiles are saved.
The OS tiles are saved in here within there own BNG folders
I then use the following code to try and build a VRT with all the tiff files that are within that directory with a CRS of BNG 22770.

  gdalbuildvrt Vector_district.vrt data/*tif

and get this message:

Can someone please help me with this bearing in mind I have little to no experience in coding or the command line?


texts:
My inputs are in this way :

col1   col2 col3
1      a     c
1      b     c
2      c     c
2      d     b
2      e     z

And I would like to obtain :

col1    col2
1       a, b
2       c, d, e

I tried doing this using the virtual layer and in the query builder using array_agg function but QGIS doesn't seem to recognize this feature.

Ideally I would like to obtain this type of results however :

col1    col2
1       a:c, b:c
2       c:c, d:b, e:z


texts: How to generate world files for raster X?

texts:
Given a list of georeferenced rasters, about a hundred, what is the simplest method to generate world files for each?

In this case the source rasters are either geotiff with embedded georeferencing or .BIL with accompanying .PRJ, but the ideal solution shouldn't care what the format is. I'm hoping to avoid duplicating data.

I have at my disposal Arcgis 10 (Arcinfo & Spatial Analyst), GDAL Utilities (gdal v1.8.0), and Quantum GIS 1.8.


texts:
I want to develop a webportal that will show the QGIS generated images on dynamic website.
For this reason I need to run gdal_translate command inside php page. but I don't know is it possible or not. If possible than please let me know how to do this and if not possible than please give any alternative solution in this regard.


texts: Create a function with many parameters and apply it over a list in R

texts:
I have a list which includes 2 data.frames and would like to apply some analysis to each data.frame. I performed my analysis to a single data.frame and I want to do the same steps over my list. I am a bit confused on how to do that but I found that the easiest way to do it is by creating a function and then to use lapply. My code for the single data.frame is below

setwd("mydir")

#import .txt file
d = read.table("vourlionas.txt",
               header = FALSE,
               sep="",
               col.names=c("x", "y"),
               fill=FALSE,
               strip.white=TRUE,
               stringsAsFactors = FALSE)

#remove the 1st row of the .txt file
d2 = d[-1,]

d3 = substr(d2$x, 1, 6)
d4 = substr(d2$y, 1, 7)

df <- cbind(d3, d4)

df1 <- as.data.frame(df)

b = apply(as.matrix(df1),2,as.character)

#rename the column names of the matrix
colnames(b) [1] = paste("x")

mode(b) = "numeric"

b2 <- as.data.frame(b)

ch <- chull(b2)
coords <- b2[c(ch, ch[1]), ]

sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID=1)))

proj4string(sp_poly) <- CRS("+init=epsg:2100")

sp_poly_df <- SpatialPolygonsDataFrame(sp_poly, data=data.frame(ID=1))

test1 = spTransform(sp_poly_df, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

writeOGR(obj=test1, dsn="mydir", layer = "r"
         driver="ESRI Shapefile", overwrite_layer = TRUE)

What I want is the above code, from the point remove the 1st row of the .txt file to be included in my function and to apply it in each element of my list.
What I have managed to do so far is to insert my elements as a list

temp = list.files(pattern="*.txt")
myfiles = lapply(temp, read.delim, sep = "")

and to write the lapply code

lapply(myfiles, testfun)

The files of the list can be downloaded from here enter link description here