texts: Installing GDAL with Java Bindings (gdal.jar) on Ubuntu 12.10

texts:
I am trying to build a postgis importer for mutliple shapefiles. After geotools took orders of magnitude longer than shp2pgsql to import my data, I decided to attempt to use the java bindings to GDAL (gdal.jar) for these purposes in the hopes that it will be much faster than geotools with an added ability to import from file geodatabases.

I built and installed GDAL using the following steps:

cd /usr/local/src/gis
wget http://download.osgeo.org/gdal/gdal-1.9.2.tar.gz
tar -xzvf gdal-1.9.2.tar.gz
cd gdal-1.9.2

./configure \
--with-python \
--with-poppler \
--with-pg \
--with-curl \
--with-geos \
--with-jpeg \
--with-png \
--with-expat \
--with-xerces \
--with-java=/usr/lib/jvm/jdk1.7.0 \
--with-jvm-lib=/usr/lib/jvm/jdk1.7.0/jre/lib/amd64/server \
--with-jvm-lib-add-rpath=yes \
--with-mdb=yes \
--with-fgdb=/usr/local/src/gis/FileGDB_API

make
sudo make install

cd /usr/local/src/gis/gdal-1.9.2/swig/java
#nano into java.opt  and change JAVA_HOME to /usr/lib/jvm/jdk1.7.0
make

I then run a test class that uses gdal.jar and add the above location to the library path so that java can find the shared objects:

-Djava.libarary.path=/usr/local/src/gis/gdal-1.9.2/swig/java

When I run the code I get the following error:

Native library load failed.
java.lang.UnsatisfiedLinkError: no gdalconstjni in java.library.path
Exception in thread "main" java.lang.UnsatisfiedLinkError: org.gdal.gdalconst.gdalconstJNI.GDT_Unknown_get()I
    at org.gdal.gdalconst.gdalconstJNI.GDT_Unknown_get(Native Method)
    at org.gdal.gdalconst.gdalconstConstants.<clinit>(gdalconstConstants.java:12)
    at elephantshape.TestGDAL.main(TestGDAL.java:22)
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:57)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:601)
    at com.intellij.rt.execution.application.AppMain.main(AppMain.java:120)

What else do I need to do to get java to recognize the shared binaries that I have built and get gdal.jar to not explode?


texts: How to handle "Error: Town not found" when updating an existing field with Field Calculator?

texts:
I want to update a field in field calculator.
update existing field is checked.
Column to be updated : Column2

CASE
        WHEN "COLUMN1"  IS NOT NULL
        THEN "COLUMN2" = "Town"
        ELSE "COLUMN2" = 0
END

Basically if Col1 has data I want Col2 to be town.
Error: Town not found.

What am I doing wrong?


texts: Creating dynamic inset map for Data Driven Pages using ArcMap

texts:
Included below  are two screenshots.

The one at the bottom is an example of an old inset map previously created which shows the extent of each data driven page.

The image  above that is a screenshot of my current predicament. Because this new project covers a much larger area and has many more pages (180 vs 35) the small inset map in the legend has to be zoomed out much farther to cover the whole area which makes it nearly impossible to see the Page Numbers when everything is actually printed out. The large purple quadrants are a viewport and each yellow box is a data driven page.

My supervisor suggested I preface the maps with a layout showing 4 different extents. Then depending on the data driven page, the extent of the inset map would just be one of these four squares and would have text denoting which section it is. This would show just a fourth of the map rather than the full extent which would allow the client to actually see the page numbers in the inset map.

I am trying to think of a rough way to approach this.

I was thinking I could add an attribute column for each of the DDP rectangles and specify which quadrant it lies within. Then I could add dynamic text in the Legend to denote which quadrant is being shown

I think my main problem lies with making the inset map actually change extent from quadrant to quadrant. It almost seems like I would need another set of data driven pages for these quadrants in the inset map, but I do not think that is possible.


texts: .net arcobjects draw polygon(rectangle) can not start at the preferred location

texts:
I am using the sample code (draw a polygon on the screen) provided by Esri but it seems I cannot start the drawing from a location I want, every time it starts at wherever the function(or tool) located on map.

Anybody know how to draw a polygon as desired? Any example there?


texts: I'm using the rasterVIS package in R to create maps with levelplot, which works fine. But I cannot figure out how to add a spatial scale bar to the maps. Does anyone know?


texts: QGIS 2.18 QuickMapServices basemap labels shrink when exported to PDF from Print Composer

texts:
For the love of all that is good on this Earth does anyone have a fix for this? When I export my print composer map in anything greater than 120 DPI the street labels from the ESRI Transportation basemap shrink to unreadable sizes. This also happens with Google Maps layers and I'm going to guess others as well.

The above image is how it looks in the print composer, but when I export as a PDF in anything over 120 dpi (In this case 300) I get the following:

sooooooo tiny.

This just doesn't make much sense to me. Does anyone have a legitimate fix? I see some similar questions around here but with no good answers.

Also, fair warning, my experience with coding is limited to a college class that I dropped after 5 days... So please, if you have an answer, don't make it too complicated. I have seen stuff about XML files and XYZ tile thingies and none of that has worked for me.... I just keep getting blank basemaps and/or errors.

To anyone who can solve my issue: you will be rewarded with my undying gratitude.


texts: How do you calculate the area of a Home Range in OpenJump - k-LoCoH

texts:
I have run k-LoCoH on OpenJump and am attempting to calculate the area of the home range. It is on a shapefile layer of the reserve the animal's GPS points are in (CSV). I have followed the manual to generate the Attributes and added Area and Perimeter to the table as explained. It generates a set of numbers for each polygon but I have no idea how this relates to area?

For example, one of the Attribute lines under the Area column has the number 8.127215840943213E-8. The largest polygon (the one I actually want the area of as it is the entire area) has the smallest number (0.013416734281735664). How can I convert this to a meaningful area estimate? Do I need to add the areas together? And how does this link to the shapefile of the reserve borders it is displayed with as they are not linked? Does this automatically calculate area as it is geo-referenced?

Apologies for all the questions but this has been driving me crazy for the last two weeks...


texts: PyQGIS: Problems with deleting attributes from Geopackage-Layer

texts:
I am struggling with a rather simple task.
While trying to delete an attribute field from a Geopackage-Layer (Polygon and Multipolygon) I encounter results I can't understand.

My Geopackage a.gpkg contains 3 layers (test_p [Polygon], test_mp and test_mp2 [both Multipolygon]).
Each has two Features and three attributes.
They all have a projection [EPSG:25833].

I am not working in the Python shell in QGIS but in Spyder 4.1.2.
Here is my code in Python 3.7 (Windows 64) for trying to delete one of the attributes in each file:

from qgis.core import *
QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.10/apps/qgis-ltr", True) # Path to qgis
qgs = QgsApplication([], False)
qgs.initQgis()

####################################################
## Delete a field
####################################################

# Path to Vector Layer in Geopackage
path_list = ('E:/a.gpkg|layername=test_p',
             'E:/a.gpkg|layername=test_mp',
             'E:/a.gpkg|layername=test_mp2')

for path in path_list:
# load Vector Layer
    vl = QgsVectorLayer(path, '', 'ogr')

    # define Data Provider
    pv = vl.dataProvider()
    # delete Field with index 2
    pv.deleteAttributes([2])
    # update fields
    vl.updateFields()

# Exit Qgis Environment
qgs.exitQgis()

Upon examining the results in the QGIS-Browser only the layer test_mp2 was altered free of error.
The other layers show in their properties no CRS and an empty extent, although the feature count is still 2.
Yet, there are no attributes.
In addition to the corrupted files two datasets of the "no geometry"-type appeared with the filename extention "_ogr_tmp". Their attribute tables are correct, as they don't show the deleted column.

It appears to be the last handled file that results in fine order. (I checked by changing the order in path_list.)

So, if someone could explain this behaviour and suggest a solution to my problem, that would be great!

Maybe I am missing some basic understanding as I just recently started to use PyQgis and am a self-taught user of python and QGis.

Best regards
Felix


texts:
I've imported a row from a Postgres DB using sqlalchemy.

One of the column types is GEOMETRY from postgis and upon import

>>>> type(row[0].geom)

returns

class 'geoalchemy2.elements.WKBElement'

I convert to a shapely polygon (the object type I need)

boundary_geom = row[0].geom
shapely_geom = to_shape(boundary_geom)

>>>> type(shapely_geom)
<class 'shapely.geometry.multipolygon.MultiPolygon'>

I'm trying to check the crs before using the object but it does not appear to have one.

I checked in the database

select ST_SRID(geom) from public.inner_london_boundary;

returns

27700

for the row in question

27700 is British national grid which is what is expected.

The Well Know Binary data passed to the function should contain the SRID.

edit: i just checked and it does:

>>>> boundary_geom.srid
27700

Why is the shapely object missing the projection metadata?

Am I just checking for it incorrectly?

do I need to set it manually no matter what?


texts: ArcGIS JavaScript API won't recognize ArcGIS online material that isn't shared publicly

texts:
I am trying to add KML Layers to the ArcGIS JS API. I have the KML layers hosted on the ArcGIS Online server. This data that I am working with contains some sensitive information and I am trying to use them within a web application that I am creating using JS (includes jQuery), CSS, and HTML. In order for me to see these layers within the web app it appears that they need to be shared publicly on ArcGIS Online. I do not want to do that because of the sensitive content of them, however I am allowed to share them within my Organization or Group on ArcGIS Online.

Unfortunately, the API doesn't recognize my permissions within these groups, and won't display the KML layers. Is there any way that I can get the ArcGIS JS API to recognize my credentials and allow me to view the layers without sharing them with everyone (public)? In other words is there any way that I can leave them shared within the Organization and/or Group and still see them within my web application using the API?

Additionally I am using Google Earth API that has the same issues, a solution to that would be helpful too.


texts: Do the recent developments in the Google Maps API replace necessity of OpenLayers?

texts:
I think the main drawback of the Google Maps API is that we can't add our own vector data or raster imagery as one of the layers in the map. Is this correct? Or is there any alternative in the Google Maps API.
If I want to add some more features to the existing features, can I do it in the Google Maps API for free?
What are the advantages of the recent developments of OpenLayers or any other client-side scripts over the Google Maps API?
If I use the Google Maps API, I think I need not require any database, webserver etc...because everything - the creation of features - is done directly using Google Maps API. Is this correct? What should I prefer to use?

I am asking all this because I want to do a project based on webGIS and I am thinking of providing a webservice for my city. Can Google Maps API alone do that or would it be better handled by OpenLayers or LeafLet, or something else?


texts: How to find out coordinate system to use pre-1940 Soviet map as raster

texts:
I am trying to use this map as a raster but I am not having any success in QGIS. I tried in Georeferencer to open raster and used these settings:

Transformation: Thin plate spine
Resampling Method: Nearest Neighbor
Target SRS: Pulkovo 1942 (and others)

Then I input five points at major lat/lng intersections and run georeferencing. But it comes out like this:

Am I doing this right? Is there a way to find out the correct coordinate system without tons of trial and error?


texts: Error in exporting attribute table into ArcGIS Desktop?

texts:
I have an attribute table of about 2.5 million rows and 10 columns that I am trying to export as dbf or text file but i keep getting the error message,

  An error occurred exporting the table, a column was specified that
  does not exist

How do i resolve this using ArcGIS Desktop 10.2?


texts: Processing feature string fields to raster

texts:
I am using the GDAL python bindings. I am trying to rasterize a layer from an arcinfo coverage file, but the values I need to rasterize I have to calculate from string fields. OGR can not write to arcinfo coverage files and, in any event, I do not want to save the newly calculated fields except to the raster. Would the following work?

Load coverage layer
Create OGR dataset of type memory
Copy layer to memory dataset
Iterate over features in memory dataset, calculating new field and adding to feature
Rasterize memory dataset

Or is there a better (or more simple) method I am overlooking?


texts:
I want to integrate GeoWebCache with GeoServer.

If a request( with a bbox say bbox-A) is cached in GWC.

Next,if a different user sends another request with a [bbox-A +- 0.00000001], will the cache be able to provide?

I believe that OpenLayers divides the screen area into a matrix (defined by image height, width) and makes a WMS request for each cell. Each request has a bbox parameter.
If the GWC is seeded with a much larger area, It internally will be saving the tiles with a different bboxes.
If there is a delta between the bbox of GWC and the request, will GWC be totally missed ?


texts: Adding own logo to the bottom of the cesium application

texts: I would like to add my own logo at the bottom of the cesium application. I would like to know which file in cesium should use to make this possible.


texts: QGIS: How to find a line with defined slope starting from a point over a DEM?

texts:
If I have a DEM as raster layer and a point with defined position in a point layer,
how could I generate a line with defined slope of e.g. -1%? The end point can be anywhere and should be defined by the line length and the constant downward slope of 1%. The line should have nodes every one or two meters, so that it follows the terrain with that accuracy.

The goal is to construct an irrigation ditch that almost follows the contour, but is defined by an optimal downward slope. If I can design it this way over a DEM with resolution <1m, I can stake it out with an RTK positioning system and have it dug accordingly.

Any ideas for a workflow or code?

I'm working in QGIS 2.18.


texts: Making st_buffer and visualizing it in QGIS?

texts:
so i have a road and i need to make a buffer of 30 meters around it. It is a road made out of few segments.
so i was thinking:

    SELECT a.ref St_buffer(a.geom,30)
FROM osm_roads a
WHERE a.ref='A306'

Am i missing something?

Do i have to do a union on the segments, because i need the road as a whole?

How do i do that?

How can i visualise it in qgis via PostGIS?


texts: How to make sure that geometries copy the relationship of corresponding tables?

texts:
A make a database in PostgreSQL (9.2.2)/PostGIS (2.0.1). I'm new to it, so I'm not sure how to make sure that geometries corresponding to main table (polygons or points) cover the geometries corresponding to the dependent table (always points).

In ArcGIS, I would simply apply a topology rule, so I tried to learn more about topologies in PostGIS, but what I found only confused me. Now I don't know how to do it through topology and I'm even not sure whether this way is a good solution, since most PostGIS topology application focus on other topology rules than I need.

Another alternative I found was a SELECT like this:

SELECT gid, the_geom
FROM foo_geometry
WHERE ST_CoveredBy(foo_geometry, bar_geometry)='f'

where foo_geometry corresponds to the dependent table and bar_geometry to the other one. I can also use a cursor for this.

What solution is the best provided that
a) as a beginner I prefer solutions easy to understand and implement
b) there won't be many geometries and the validation will be done just once in a time as a maintenance process, so performance is not crucial here
?


texts: Geoserver cascading and fetching projection

texts:
I tried to use cascade of WMS from this portal in order to cache wms layers

http://geoportal.dgu.hr/wms?request=GetCapabilities

Native for this portal is EPSG:3765, but EPSG:31275 is also supported
But while direct access to WMS returns correct result:

http://geoportal.dgu.hr/wms?layers=DOF&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A31275&BBOX=5578043.9690343,5053103.9883878,5585211.9651635,5060271.9845171&WIDTH=256&HEIGHT=256

cascaded WMS  with exactly the same parameters returns all-white image.
And in documentation of Geoserver it says that while cascading if the requested projection exists on remote server it will serve it from remote server, otherwise it will reproject it.

What could be going wrong, why isn't geoserver just returning correct image from remote service?

Geoserver version used is 2.2.


texts: How can I dock ArcToolbox in ArcMap 10 custom template and make it stay?

texts:
Here is the situation:
Because UI Controls cannot be saved to Normal.mxt like they used to I had to create an empty document with all my custom VBA code that is mostly called from the local document-level UI controls. (Until I have time to get proficient with .NET things will have to run this way)
All is good although it would have been nice to have Normal.mxt-level UI Controls back...
Whenever I need custom functionality i start creating the new document within my Custom Template. One annoying problem I have is that I have to dock ArcToolbox every time i open a document created from the custom template. It will just not stay docked whether I save the custom template with it docked or the Normal.mxt. Normal template documents are fine, ArcToolbox is where i want it, just the custom template is misbehaving.
Any ideas?
Thanks


texts:
My .NET/OpenLayers/GeoServer/PostGres-PostGIS application allows users to create custom queries to get data from my PostGIS tables.

The queries are called in .NET (not javascript) using a web method  using the following SQL query:

SELECT "tract_name", ST_AsGeoJSON("geom") FROM a_postgis_table;

The query works fine when I run it in the PostGres SQL window.

How can I pass the  PostGIS results, including the GeoJSON column, back to javascript/openlayers from the .NET web method so that OpenLayers can render the vector layer in the map? There has to something really simple I am missing here. Seems like this should be an extremely common thing that people need to do.


texts: Dynamic text colour in ArcGIS Pro Table Frame

texts:
Has anyone figured out a way to alter the text properties (colour, boldness, font, etc.) of a table element based on a specific property in ArcGIS Pro?
For example, I want to add tables to a layout in Pro and have certain cell text highlighted if they are above or below a certain value or colour a specific way if the value in the cell exceeds a threshold.
I'm using the most recent version of ArcGIS Pro where table elements have been added and can be dynamically linked to the extent of the map.


texts:
I want to export the SRTM 30 m elevation data and slope for the country of Indonesia. Here is my script using GEE.

However, on exporting the tiles to my Google Drive, I get this error message- A region must be a GeoJSON Polygon or LinearRing. Got: 'MultiPolygon'.

How do I deal with this i.e. how do I export the elevation (and slope) tiles for Indonesia to my Drive and then download it to my computer?

// Load the SRTM image.
var srtm = ee.Image('USGS/SRTMGL1_003'); //Calling SRTm 30m images

// Apply slope algorithm to an image.
var slope = ee.Terrain.slope(srtm);

//Scale and Projection
var scale = srtm.projection().nominalScale();
print('SRTM scale in meters', scale); //see Console tab in the right section of screen- 30.92208 m
print('Projection, crs, and crs_transform:', srtm.projection()); //EPSG:4326 is WGS 1984, look at right section of screen
// in the Console tab

// Load countries border
var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw');
var justcountry = countries.filterMetadata('Country', 'equals', 'Indonesia'); //change country name here
Map.addLayer(justcountry)

//Displaying slope for country
var elevcountry= srtm.clip (justcountry);
var slopecountry= ee.Terrain.slope(elevcountry);
//var slopecountry_radians= slopecountry.divide(180).multiply(Math.PI) //converting to radians 'divide the slope by 180, multiply the result of that by π
Map.addLayer(slopecountry); //to view

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: elevcountry,
  region: justcountry.geometry().bounds(),
  description: 'elev',
  scale: 30, //pixel size specified
  maxPixels: 1e12
});


texts: How legitimate is smoothing raster data that has been converted from polygons?

texts:
I am interested in applying this methodology to census block data after conversion from a raster:

http://rutgerscps.weebly.com/uploads/2/7/3/7/27370595/censusdatasmoothing_brief.pdf

Basically, the steps are:

Calculate the average distance between polygon centroids - Average Nearest Neighbor tool in Spatial Analyst
Convert Features to Raster
Smooth the outputting raster data by applying the focal statistics tool using the radius calculated in the Average Nearest Neighbor tool.

It seems that, although this may make sense intuitively, because sharp transitions in population count or housing units don't always occur, they may occur sometimes.  Though the result is 'nicer to look at', in that it's more a heat map than a rigid rasterized layer, is this workflow statistically defensible?  What tests (other than the simple one provided in the article) could you run to test out the smoothed result?


texts: Understanding the cell size argument in the kriging tool, ArcGIS

texts:
I'm not quite understanding the cell size values that a given raster has. I've been doing some kriging interpolation on very dense point data and if I leave the cell size parameter blank, it's cell size is like an extremely small number, like e-006.

What does this extremely small cell size number mean? I always known that the cell size of a raster is the size of the cell that each value of the raster represents spatially. But with such a tiny number, doesn't seem to make sense.  When running the Kriging tool and if I enter a whole number, like 1 for the cell size, it just errors out. If I zoom into the raster and measure the size of the cell pixels, they look like there around 3 ft by 3 ft in size. Am I confusing cell size with pixel size here? If so, what determines the size of the pixels?


texts: Import more than one GeoJSON into Leaflet map?

texts:
I've searched far and wide and tried a number of different options but I cannot get more than one GeoJSON into my Leaflet web app. I can get one, and once I try anything else, the map itself will go blank!

Below is the code I have for the input of the first GeoJSON.

Below is the line I tried to add, I've tried 5 different solutions with no success, so I would be grateful for any help at all.

    $.getJSON("turk.geojson", function(data) {

    var geojson = L.geoJson(data, {
      onEachFeature: function (feature, layer) {
        layer.bindPopup(feature.properties.Area_Name);
    }
    });

    var myStyle = {
    "color": "#ff7800",
    "weight": 5,
    "opacity": 0.65
};

    var geojsonLayer = new L.GeoJSON.AJAX("data/tur_route_geo.geojson" , {style:myStyle});

    var map = L.map('my-map')
    .fitBounds(geojson.getBounds());
//    .setView([0.0,100], 2);

    basemap.addTo(map);
    geojson.addTo(map);
    geojsonLayer.addTo(map);
  });
};
</script>

#error
127.0.0.1 - - [10/Jul/2019 10:03:39] "GET /tur_route_geo.geojson HTTP/1.1" 304 -
127.0.0.1 - - [10/Jul/2019 10:03:52] "GET / HTTP/1.1" 200 -


texts:
Currently, I would select a bunch of polylines in a fgdb, using Selection > Selection by Location. I will then be using a "AOI" (Area of Interest) polygon to help me select the polylines in a fgdb that are within AOI (usually, i will choose "are completely within" option), and these selected polylines will appear on one of my mxds.

I would like to have a script that can help me do this: running the selection of polylines using AOI polygon as a reference. So, whenever i open up that mxd, i will be able to see the updated polylines on it.

FYI, the feature class in fgdb will be updated monthly (names of fgdb and its dataset remain unchanged), and there might be a increase or decrease in the number of polylines every month (there are to be captured as long as there are within the AOI.


texts: What causes the GPS offset/shift in China?

texts:
I've spent close to two full days trying to understand the nature of the infamous GPS shift in China. After sifting through the noise, two recurring symptoms are that:

GPS coordinates (WGS-84) plot poorly on Chinese maps, with an offset ranging from 50 to 1000 meters.
Google and Apple Maps in China have severe misalignment between their street maps and satellite imagery (live example)

Various sources make contradictory claims:

Cold War-era regulations mandated introducing random offsets into GPS chips due to security concerns
protectionist measures aimed to favor GPS devices manufactured in China
China simply uses a different projection system, and the offsets are not random
GPS chips manufactured in China offset the coordinates they return, in order to match approved Chinese maps, which must follow the GCJ-02 datum

What is the actual story? This looks like the type of pervasive problem that merits a Wikipedia page, but I haven't been able to find anything authoritative yet.

The two questions suggested as possible duplicates don't answer exactly what the cause is; rather, they describe symptoms.


texts: What datastructure is a good choice for objects in a map-editor?

texts:
I am planning to implement a map-editor. But I don't really know what data structure are a good choice for storing the map-object in.  When the user pans (scrolls) the map, there will be much looking up of what objects are in the visible area, usually by having an min_x, min_y and max_x, max_y coordinate pair.

I have read that R-trees are good for this type of use, but they seem to be complicated to implement.

So far I have just used a linked list that I have to loop through very often; it's not optimal, but it works. Are there any better choice of data structures?


texts: Validity of ordinary kriging of polygon centroids

texts: I'm in a discussion with a colleague who believes that kriging should not be used to interpolate a surface from polygon based data (i.e. with kriging applied to centroids), because this effectively translates the spatial uncertainty inherent in any polygon centroid into data error in the kriged surface.  I'm sure he's right in principal but cannot help thinking that the uncertainty-error question should not be an absolute argument against the method give that uncertainty is (a) unavoidable and (b) a function of zone size and the degrees to which the data has already been aggregated.  I'd be interested to hear views from others, particularly anyone with good experienced of kriging interpolation.


texts:
I have followed the instructions on this blog post on adding streetview images to the infowindow:
http://blog.cartodb.com/a-picture-is-worth-1-000-attributes-infowindows-with/

I can do it without a problem and the images work well, however while the SQL query is running I cannot add new points to my map while in the CartoDB editor.

Is there a way around this. I would like to show the images but also be able to add points.

I am using the SQL query below. Would the solution be to change it to a query which actually writes the data into my table permanently for every new point I create, rather than just querying it? If so what would that script look like?

(Sorry about my limited SQL knowledge)

SELECT *, 'http://maps.googleapis.com/maps/api/streetview?size=300x190&location='||
ST_Y((the_geom))||','||ST_X((the_geom))||'&sensor=false&fov=110' as image
FROM table_2012_entire_tree_down_complaints_sandy

If this is possible then would a Google reverse geocode also be possible, utilising a similar script?


texts: Is there any simple Map Labeling Algorithm?

texts:
I was looking into how to Label a few features on a small Map, and want to label them on the fly.

I would be adding a few features (a maximum of around 6 to 7) as graphics on a web-Map which has only a hypsometric rendered DEM. I would like to label these graphics and make sure that the labels do not overlap

I have looked at quite a few papers, but the algorithms look quite complicated. Is there any simple Map Labeling Algorithm that gives reasonable results?


texts: How do I generate thumbnails images from lots of shapefiles?

texts: A colleague and I recently released an Open Source project called Mapoteca which is a set of Python Scripts that integrates with a GDAL/FWTools installation allowing a full walk-through of directories looking for spatial data and extracting extent, geometry, number of elements, projection and user writing it into a csv file. We are working on the part that describes the DBF files associated and we wanted to generate a thumbnail image of each shapefile, do you know of a way to generate these thumbnails in jpeg or png using Python?


texts: How to combine only overlapping parts of multiple shapefile layers?

texts: I have 5 shapefiles. 4 of roads and one of 3 adjoining cities. How can I combine somehow the areas where the street and the cities shapefiles overlap, omitting streets not lying in those cities? The colored areas in the lower left are the cities.


texts: get date from image using python API in Google Earth Engine

texts:
I am trying to extract the date as a string from an image in google earth engine, but can't seem to get it right. This is the code I have so far.

date = ee.Date(image.get('system:time_start'))

This returns an ee.date object in Python, but I can't find the right way to access the object. Date.getInfo() returns an error.

Does anyone know how to do this?


texts: Replicate Import Symbology Matching Dialog in ArcObjects?

texts: In Arcmap I can render a layer with same symbolgy applied to different fields by way of re-importing the .lyr file and choosing a Value Field.  How can I replicate this  functionality in ArcObjects?


texts: Improve speed of PostGIS query for counting locations

texts:
I have one table that imported from shp file (projection - wgs84- polygon).

The table (locations_list)has 5 fields: id, geom, pin, shape_star, shape_stle, isEmpty.

I created spatial index on the table with a geometry column (geom)

pin is the location id and it is unique.

I have more than 15k locations(pin) and each location either empty or not(isEmpty takes values 0 or 1)
I want to count for each location(pin) in the table how many empty location around it within 500 feet.
I used the following query to  do that:

Create table as count
Select a.pin
sum(a.isEmpty) as count
from locations_list a
join locations_list b on a.pin != b.pin
and
ST_Dwithin(a.geom::geography,b.geom::geography,152)
group by a.pin

This query is under testing and it takes hours to complete.
Is there any advice to improve the performance?


texts: What can I do with this unzipped GDB file in QGIS?

texts:
I have a zip file that I got from the City department of information technology. Unzipped, it is a directory called DOITT_ZIP_CODE_01_29JUL2009.gdb containing a bunch of .gdbtable and .gdbtablex files. I'm not super seasoned in GIS but I'm pretty sure I ought to be able to open this in QGIS. I can't figure out a combination of parameters that will let me open it, however.

I really just need to get the data out, so I'd be happy to use something like GDB Flee but it can't seem to read the .zip and the .gdb is a directory.

Note: These instructions (Migrating geodatabase data into PostGIS without Esri apps?)  suggest I could use a command line tool, ogr2ogr but mine doesn't seem to have FileGDB support.

What is my next step?


texts: Is there a way to get Openlayers' KML reader to deal with tessellation?

texts:
I'm trying to use OpenLayers' KML reader to display a KML file on a map, but the KML itself has LinearRings that have a tessellate value of 1.  I realize that this is used by google earth to snap the object to the ground for a more accurate drawing, but it causes the KML to be useless as far as the OpenLayers KML reader is concerned.

So my question then is, what can be done to convert these tessellated regions?  I'm assuming there's not simply a setting I can change in OpenLayers to have it properly display the regions without tessellation.


texts: Solutions for fast operations on large GIS datasets

texts:
Possible Duplicate:
  Acquiring ArcGIS-like speed in Postgis

I have been using PostGIS 2.0 for all of my GIS operations. PostGIS is an excellent GIS database, but no matter how I format my queries or design my tables I still have queries that take many hours or days to run. These same operations, when carried out in ArcGIS would take orders of magnitude less time to run comparatively. It is very depressing and stressful to just sit there and wait hours for simple queries to complete when you know they could be completed in much less time and with fewer chances for error. ArcGIS itself is not an option because I don't just have 75,000 dollars to spare for a windows-only program with strict restrictions and an ancient version of python.

Consider a simple point in polygon operation on a parcel dataset of a few hundred thousand parcels with a few hundred thousand points. This can take 45 minutes to an hour (with everything correctly indexed and using st_intersects); whereas in ArcGIS it would take less than a minute.

I have been told that this is because postgis is meant to handle operations on individual features rather than layers of features. In my experience though, QGIS does not perform operations much faster than in postgis and is much more subject to random crashes and generally seems unstable. The other thing that bothers me is that I have 3 local machines, each with an i7 processor and 8-16 GB of ram just sitting around doing nothing when they really are needed for helping out with GIS processing --but nothing I have at the moment can utilize them and I seriously doubt for that matter that my system even uses my single machine efficiently for these purposes.

I am very interested in clusters and distributed map reduce for GIS with NoSQL --from what I can tell it seems like these technologies are probably what I need to build anything that is significantly faster than what is out there. The problem is I really have almost no experience in these areas and am a very new programmer with only a working knowledge of python sql and a little java. Furthermore, the only GIS programs I have ever used are postgis, QGIS and ArcGIS.

The various options seem infinite so I am looking to the GIS community for some direction on where a beginner with no networking experience like myself should start and what technologies I should use to build a faster --or improve an existing-- GIS platform.

My requirements are as follows:

Significantly faster than postgis for raw data processing on large datasets.
Is supported in that it uses standard modules like GEOS or JTS, mongoDB etc. --I don't really want to build complex mathematical/geometric algorithms on my own as I lack the background to do this (at present, although I am slowly making an effort to get better at such things). In other words I want to use software that isn't extremely obscure and has a vibrant user/dev community.
Works on Ubuntu server 12.04+ (I suppose debian would be fine too).
Uses my unused machines on a local network. --This is not strictly a requirement as ArcGIS doesn't use multiple machines to get blazing speed; however, since they are just sitting there I might as well use them for something if at all possible.

I anticipate that with any of these options I will have the following two questions based on my research:

How do I store geometries in a standard way that I can easily convert to postgis or shapefiles.
How do I create an efficient index on the data.

Thanks.

EDIT:
Due to all of the follow up questions in the comments being about things I could have possibly done wrong in postGIS, I have posted a question that more specially addresses the problems I have with postgis, so that this question can specially refer to what it would take to create a better platform


texts: Looking for Fortran code calculating drainage basin area

texts:
I am working in a project of glacier mass balance.
For my project I want to integrate the area of the drainage basin in a glaciological system model (GSM). As the drainage area changes with time due to topographic evolution (ice accumulation or melt), I need to be able to integrate some code which calculates this area in the GSM. Would anybody know where I could find some code to compute the drainage area of a region?


texts: Creating geodesic hexagonal global grid in QGIS?

texts:
I'm looking to make a geodesic hexagonal global grid based on a icosahedron. Such a grid would need to have 12 pentagons in it as well to be able to fit a sphere.

I know mmqgis can generate grids, but these grids are plain flat grids, not geodesic ones, and they do not map to a sphere without extreme distortion, which I'm trying to avoid.

What I'm trying to make is something like this:

Ideally, this would be done on windows, but Linux is also an option if there are no tools available for the former


texts:
I'm working in Qgis with python console.
I have a file .asc whose structure is:

Ncols 49

Nrows 52

Xllcorner 2611187.75

Yllcorner 4336914.199

Cellsize 1000

NODATA_value -9999

0.5081 0.679 -9999 ...

I have to find the Max and the Min value in the file and I tried:

GridFile = open('filepath \ file.asc', 'r')
MaxValue = max (GridFile)
MinValue = min (GridFile)

But the procedure returns the following error

ValueError: max () is an empty sequence
ValueError: min () is an empty sequence

What I understand is that I have to skip the first 6 lines and find the maximum and the minimum of successive values, but I do not know how to do it.


texts: Multiple-selecting using Shift key in Project pane of ArcGIS Pro?

texts:
Using ArcGIS Pro 1.2, I just tried to use and the Shift key to multiple-select a set of shapefiles in a folder connection within the Project pane.  I thought I had done this before without a problem, but I am now thinking that I must have imagined that because if I left-click on one, and then hold the Shift key to left-click on another only the second one is highlighted.

I can use the Ctrl key to multiple-select, which is fine for a few items, but not when I may have many.

If the Shift key is expected to be able to multiple-select in the Project pane, is there a known cause for why I might not be able to?


texts: QGIS time manager plugin - animated trajectory: how to keep focus on a moving object

texts:
I have just animated 4 different tracks from boats on the Ionian and Adriatic Seas by using time manager plugin. Everything works fine now (since I swiched over to shapefiles from GeoJSON).

The area covered by the tracks spans about 700km and I would like to render a video where time and speed labels are visible. So I would like to define a certain scale and keep one boat centered within the display during the time manager loop. Is there an easy way to do that?


texts:
The problem I'm trying to solve is finding layout of street lights on a street. The street lights can be arranged in the following layout: one-side, two-side opposite, staggered, or median, etc. as shown in the picture.

I have latitude, longitude of street light points of several cities and I want to build a classifier that can determine the road lighting layout.

I thought of the following approach so far but I'm stuck at what should I do after a step.

Reverse geocode all the street light points using Nominatim and assign each one to a street_id using street name as well as get the OSM way id.
Use the OSM way id to obtain the polyline of the street.
Find out which points lie on which side of the polyline using Linear Referencing Toolbox in ArcGIS.
Find the distance between poles by calculating average distance to nearest neighbor on each side of the road.

I'm stuck here.


texts: Errors in creating spatial index on SDO_GEOMETRY column in Oracle 11g

texts:
I am using Oracle 11g (Oracle Database 11g Enterprise Edition Release 11.1.0.7.0 - 64bit Production) and I have a table that has column of type SDO_GEOMETRY. I want to create a spatial index on this column (NEW_SHAPE). My query is as following:

CREATE INDEX GIS_GEOM_SRID3857_SPTIDX ON GIS_GEOM_SRID3857_LOOKUP ( NEW_SHAPE )
INDEXTYPE IS MDSYS.SPATIAL_INDEX;
/

and I am getting the following errors:

Error starting at line : 23 in command -
CREATE INDEX GIS_GEOM_SRID3857_SPTIDX ON GIS_GEOM_SRID3857_LOOKUP ( NEW_SHAPE )
INDEXTYPE IS MDSYS.SPATIAL_INDEX
Error report -
ORA-29855: error occurred in the execution of ODCIINDEXCREATE routine
ORA-13249: SRID 3857 does not exist in MDSYS.CS_SRS table
ORA-29400: data cartridge error
Error - OCI_NODATA
ORA-06512: at "MDSYS.SDO_INDEX_METHOD_10I", line 10
29855. 00000 -  "error occurred in the execution of ODCIINDEXCREATE routine"
*Cause:    Failed to successfully execute the ODCIIndexCreate routine.
*Action:   Check to see if the routine has been coded correctly.

I have searched on these errors but have not found anything useful. What is the problem here?


texts: Qgis 2.12 Raster Calcuator gives nan raster image

texts: With the 'r.los' command I've made some tiff raster images from a dem tiff. When I open the raster calculator to do math with the new rasters, the result is always a "nan" raster. From its options I load min/max values and it gives me 1.79769e+308 and -1.79769e+308. I reclassified the 'los' rasters with integer values but no change. I used the clipper to fill no data values with 0, but nothing changed. I saw here a post which said to do the math operations like this "float(los1+los2)" but nothing changed. The "fill no data" command gives me values from 1-255. The "r.fillnuls" command corrupts the results probably because i don't know what value to put in the "Spline" parameters. Do you have any option?


texts: ESRI ArcGIS collapses geometry arrays for nearby points

texts: When calling ESRI ArcGIS Solve Route REST API (https://developers.arcgis.com/rest/services-reference/solve-route.htm), if the input stops are nearby (for example -121.31678799999997, 44.024490000000071 and -121.316788, 44.023907, the returned output considers them as a single point in the features -> geometry JSON element. If there are total n input stops, there should ideally be n-1 arrays in the output (each for a path between two stops) but if a pair of points like the above is in the stops, they 'collapse' and the output contains only n-2 arrays. How to prevent them from 'collapsing'? I have tried using the parameter outputGeometryPrecision with value 0 and then again with value 1 and value esriCentimeters for outputGeometryPrecisionUnits to no avail.


texts: How do I convert a latitude/longitude coordinate to a NOAA GFS grid cell location?

texts:
Looking at this GFS http://nomads.ncep.noaa.gov:9090/dods/gfs_0p50/gfs20170707/gfs_0p50_12z.ascii?tmax2m[0:2][124:128][279:283]

The range defined by the second and third set of brackets specify the grid coordinates, however I'm not sure how to convert between what NOAA is using and the common geographic coordinate system.

What would the NOAA equivalent be for this coordinate pair: (36.25, -78.75)?


texts: What is Login to Esri FTP site for Arc Hydro Tools?

texts:
Esri used to have a forum sticky that contained the login instructions to the FTP site for the Arc Hydro tools. They have since deleted all of their forums.

Does anyone know the login instructions for the Esri Arc Hydro FTP site?


texts: How useful is the ESRI Production Mapping extension for cartography?

texts:
I would like to find out how useful is the $9000? ESRI Production Mapping extension.  According to ESRI is helps "Standardize and centralize detailed cartographic production with tools for creating and maintaining derived data, symbology, page elements, and maps."

Can someone provide opinions on contour creation and labelling, layout features, styles (apparently stored within a geodatabase?) Area there any enhancements to the legend creation process?


texts: Restricting road access by vehicle type in ArcGIS Network Analyst

texts:
Is there a way to restrict road access by vehicle type? Some roads restrict trucks, some roads have height/weight restrictions, and some roads are restricted by the type of license plate (e.g., commercial, personal, etc.).

I am familiar with restricting road access by height and weight of the vehicles, but how would I add attributes in Network Dataset to restrict road access by vehicle type (e.g., trucks) and type of license plate (e.g., commercial)?

Is it possible to achieve with Network Analyst extension?


texts: how to create RTree index on LineString geometry

texts:
I'm working with NTS(NetTopologySuite), then I get some geometries with type of LineString, and I want to create the IntervalRTree Index on them.

Here are some codes:

SortedPackedIntervalRTree<LineString> tree = new SortedPackedIntervalRTree<LineString>();

and the Insert Method looks like:

public void Insert(double min, double max, T item)

so, how can I get the "min" and "max" Value of LineString?

also, the Query Method of the tree looks like

public void Query(double min, double max, IItemVisitor<T> visitor)

and how to query the tree index?


texts:
I want to create a Desktop project on Java which shows radiation dispersion. How can i Link java with QGIS?

Is there any tutorial site i can get reference from?


texts: ogr2ogr FAILURE: Couldn't fetch requested layer 'ùnlt'!

texts:
I am getting the following error message with ogr2ogr while try to convert into MySQL.  I have tried PGeo,Shapefile and FileGDB.

FAILURE: Couldn't fetch requested layer 'ùnlt'!

My environment is Window10 running under Paralles11 on Mac OS X El Capitan.
I have tried in FWtools and latest OSGeo4W Shell.  Has anyone seen this error before and know how to resolve it?


texts: InsertCursor fails when trying to insert XY point data

texts:
I have a text file with X and Y coordinates, as well as a third piece of information. I'm trying to create a point feature class from this text file by use of an insert cursor. I'm working with "SHAPE@XY" token, but can't figure out the exact syntax. What am I doing wrong?

Pertinent code:

#Open text file
textFile = open (inFile, "r")

#Get X, Y, and PGA
lines = [r.split (" ")[0:3] for r in textFile]

#Spatial reference
SR = SpatialReference ("WGS 1984")

#Create point feature class
CreateFeatureclass_management (outPath, outName, "POINT", spatial_reference = SR)

#Add PGA field
AddField_management (outFc, "PGA", "FLOAT")

#Create insert cursor object
cursor = InsertCursor (outFc, ["SHAPE@XY", "PGA"])

for X, Y, PGA in lines:
    #XY tuple
    XY = (float(X), float(Y))

    #Create row tuple
    row = (XY, float(PGA))

    print row
    #Insert row
    cursor.insertRow (row)

The first line prints as I would expect:

((-123.224, 38.4105), 0.02)

But the next line throws the following error:

Traceback (most recent call last):
  File "C:/Users/e1b8/Desktop/E1B8/DASH_GIS_Integration/Calculate_PGA_per_plat_Point.py", line 277, in textToFeatureClass
    cursor.insertRow (row)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 106, in insertRow
    return convertArcObjectToPythonObject(self._arc_object.InsertRow(*gp_fixargs(args)))
RuntimeError: Error in InsertRow

Three sample lines from the text file:

-123.2240 38.4105 0.02 0.01 1
-123.2157 38.4105 0.02 0.01 1
-123.2073 38.4105 0.02 0.01 1


texts: Convert 32 bit float raster file to JP2000?

texts:
Currently, i have a 32 bit float raster (.TIFF) and would like to convert it to JP2000 format. Is there a way to do this?

I have tried Copy Raster tool and Raster to Other format tool but both produce the error:

  000446: Output file format with specified pixel type or number of
  bands or colormap is not supported.

My raster is single band without color map and 32 bit float. Here is the statistic of my raster:

Minimum: 493.33380126953
Maximum: 493.33380126953
Mean: 662.78674316406
Std Dev: 49.161709494491

And my desire precision is 3


texts: Query one overpass api process multiple times

texts:
I using a local instance of the overpass api, because I need to make a lot of queries (over 100k). The solution I'm using right now is to call bin/osm3s_query for every single query.

This is very slow, because osm3s_query has to start, answer the query and shut down for every single query.

Is there a way of having some kind of daemon waiting for my queries?

I guess if osm3s_query doesn't have to start and shutdown for every single query, this should accelerate things a lot.
The "daemon" for running a web server seems to do exactly what I'm doing right now (firing up one osm3s_query process for every query). Therefore I couldn't find anything helpful in the documentation


texts: Getting coordinates of ways with Overpass API osm script

texts:
I realize there are questions similar to this one but after looking at them, I still can't fix my problem.

I would like to get all the node coordinates for the nodes along all the road ways in some area. Right now I can only get which nodes are contained in the way (w/o their coords). From answers to other questions here I think that I need to change the print mode to geom but I can't get it to work.

Here's my query:

"""
<osm-script output="json" timeout="25">
  <!-- gather results -->
    <query type="way">
      <has-kv k="highway"/>
      <has-kv k="highway" modv="not" v="footway"/>
      <has-kv k="highway" modv="not" v="pedestrian"/>
      <has-kv k="-highway" modv="not" v="path"/>
      <bbox-query s=" """ + str(s) + '" w="' + str(w) + '" n="' + str(n) + '" e="' + str(e) + '"/>' + """
    </query>
  <!-- print results -->
  <print mode="skeleton"/>
  <recurse type="down"/>
  <print mode="geom" order="quadtile"/>
</osm-script>
"""

(Note: vars s, w, n, and e are the max/min long and lat coords defined elsewhere)

When trying to convert the output to Json:

response = requests.get(overpass_url, params={'data': overpass_query}, proxies=proxies)
data = response.json()

I get the following error:

---> 37 data = response.json()
JSONDecodeError: Expecting value: line 1 column 1 (char 0)

I can provide the original query that results in a Json w/ the nodes in each way but w/o their coordinates if that is necessary, but I thought it might be overkill.


texts: Re-projecting Multiple Shapefiles using ArcPy?

texts:
I decided to re-project shapefiles from one coordinate system to another. I wrote the following script but am getting the execution Error 999999.

Here is the script:

import arcpy

# input data is in NAD 1983 UTM Zone 11N coordinate system
Input_Elk_Feature = "C:\Temp\Data\Originals\Elk\ds945.shp"
Input_ElkHuntZones_Feature = "C:\Temp\Data\Originals\Elk\ds786.shp"

# output data
Output_Elk_Feature = "C:\Temp\Projected_Data"
Output_ElkHuntZones_Feature = "C:\Temp\Projected_Data"

# create a spatial reference object for the output coordinate system
Out_coordinate_system = arcpy.SpatialReference('NAD 1983 UTM Zone 10 North')

# run the tool
arcpy.Project_management(input_features, output_feature_class, out_coordinate_system) #reprojects the data into the specified coordinate system.

TestList = ["C:\Temp\Data\Originals\Elk\ds945.shp", "C:\Temp\Data\Originals\Elk\ds786.shp"]

x=0

while x>=0:
    Input_Elk_Feature = TestList[0]
    Output_Elk_Feature = "C:\Temp\Projected_Data"+0+".shp"
    arcpy.Project_management(input_features, output_feature_class, out_coordinate_system)

    x=x+1

    Input_ElkHuntZones_Feature = TestList[1]
    Output_ElkHuntZones_Feature = "C:\Temp\Projected_Data"+1+".shp"


texts: Printing WMS map with chosen scale and WMS zoom level using QGIS?

texts:
I am trying to print a WMS  map in QGIS.

If, while in the print composer, I chose a "25000" scale, the composer chooses a WMS layer that actually shows a 1:50000 map. To get down deep at the 1:25000 WMS map in the print composer preview, I have to trick it by telling the composer to use a smaller scale, but then I don't get the scale I want.

How does QGIS choose its WMS zoom level?

I need to chose a scale AND a specific WMS zoom level. I tried the "Lock layers for map item" check box but it is not what I hoped for.


texts:
I'm trying to use pyqgis to retrieve the name of the group in which a layer with a specified name is located. My function code is as follows:

def runA4(self,inputRasterName,outputFilename,outputLayerName):
    root = QgsProject.instance().layerTreeRoot()
    print "root:"
    print root

    demRaster = root.findLayer(inputRasterName)
    print "demRaster:"
    print demRaster

    parentGroup = demRaster.parent()
    print "parentGroup:"
    print parentGroup

This gives the following output in the Python console:

root:
<qgis._core.QgsLayerTreeGroup object at 0x0000000025287730>
demRaster:
None

And the following in the "Python error" log messages panel:

2017-10-17T16:42:40 1   Traceback (most recent call last):
    File "C:/Users/Gareth/.qgis2/python/plugins\myPlugin\myPlugin.py", line 1749, in runA4D
        self.runA4(inputRasterName,outputFilename,outputLayerName)
    File "C:/Users/Gareth/.qgis2/python/plugins\myPlugin\myPlugin.py", line 1761, in runA4
        parentGroup = demRaster.parent()
    AttributeError: 'NoneType' object has no attribute 'parent'

So there's something going wrong when I try to grab demRaster using inputRasterName, which is a string containing the name of the layer as listed in the layer tree. What am I doing wrong here?

I can't help but think there must be a better way to do this - it seems like it should be a very simple operation!


texts: Locating 0-points multipart polygons using arcpy

texts:
By Running this piece of code:

import arcpy
arcpy.env.workspace = "C:\working\GEODATABASES\\tests.gdb"
fcList = arcpy.ListFeatureClasses()
for fc in fcList:
    with arcpy.da.SearchCursor(fc, ("OID@", "SHAPE@")) as cursor:
    for row in cursor:
        print("Feature {0}".format(row[0]))
        print "   " + str(row[1].partCount) + " parts for THIS PARTICULAR FEATURE"
        for part in row[1]:
            poly = arcpy.Polygon(part)
            print("   Points {0}".format(poly.pointCount))

In theory I am able to list all polygons for a particular feature class and list whether they are multi-part polygons or not. In addition, it shows the number of points per feature part.

HOWEVER, it does not do the job properly as I am dealing with a complicate feature class that actually have more feature parts than the script can identify. Even in some of them, it does not obtain the right number of points per part either.

By using the OLD reading cursor system in arcpy (not the above one), I get in some polygons a different number of points, but still there are some of them missing, therefore the script cannot access to them.

Why is that this script stops counting either points or parts sometimes when entering a multipart-polygon within my feature class?

Is it my code not well tuned for complicated multi-part polygons?

Thanks,


texts:
I'm working on a project where I want to deal with one-to-many relationships using a geopackage. I can create layers and non-spatial tables but I'm having trouble with the spatial view. I want to have many condition audits conducted on a footpath and want a view to display only the most recent condition audit but QGIS doesn't seem to be working with a spatial view in the geopackage.

Software versions:
QGIS 2.18.3,
'DB Browser for SQLite' version 3.9.1

workflow:

Using QGIS, create a polygon layer using 'Layer>Create Layer>New Geopackage Layer'
Then create a non-spatial layer using the same method and selecting the same geopackage but select the ‘Geometry type’ as 'non spatial'
Using DB Browser enable foreign keys checkbox in pragams and add a foreign key constraint to non-spatial table (FootpathCondition)
Then create a spatial view:

CREATE VIEW FootpathConditionRatings as Select fp2.fid as fid, fpc2.fid as AuditID, fpc2.AuditData as DateAudited, fp2.geometry as geometry from Footpath as fp2 INNER JOIN FootpathCondition as fpc2 ON fpc2.FK_Footpath = fp2.fid where AuditID in( select AID from( Select fp.fid as fid, fpc.fid as AID, MAX(fpc.AuditData) as DateAudited from Footpath as fp INNER JOIN FootpathCondition as fpc ON fpc.FK_Footpath = fp.fid group by fp.fid ) )

add values to gpkg_contents and gpkg_geometry_columns to make the view a spatial layer
using QGIS populate tables adding 2 polygons and 8 records in non-spatial table (representing 4 transactions against each polygon)
Using DB Browser check 2 tables' and the view's data. View returns 2 rows with the most recent date, as expected.

fid AuditID DateAudited

"2" "3" "2017-02-14"

"1" "8" "2016-05-14"

Add view to QGIS using DB Manager and check data

Outcome:
QGIS displays the view's geometry correctly however the 2 records in the attribute table are the same values. If I copy the 2 records to the clipboard and paste then these are the values I get:

    wkt_geom    fid AuditID DateAudited

    NULL

If I select a row in the attribute table both rows are selected and both geometries are also selected.

I'm wondering if anyone else is getting this behaviour with spatial views in a geopackage?


texts: Trouble displaying GeoJSON file in OpenLayers

texts:
I am noobie to OpenLayers and can not for the life of me figure out how to display a GeoJSON file. This seems like it should be easy but I have been banging my head for a day. I'm sure I am missing something obvious. I have my GeoJSON file (countries.json) on a server and basically followed the instructions in "OpenLayers:beginners guide" word for word so I don't know why the thing won't display. Here is my code:

function init(){
    map = new OpenLayers.Map('map_element', {});
    var vector_strategies = [new OpenLayers.Strategy.Fixed()];
    var vector_format = new OpenLayers.Format.GeoJSON({
        externalProjection: new OpenLayers.Projection("EPSG:4326"),
        internalProject: new OpenLayers.Projection("EPSG:900913")
    });
    var vector_protocol = new OpenLayers.Protocol.HTTP({
        url: 'ne_50m_admin_WM.json',
        format: vector_format
    });

    vector_layer = new OpenLayers.Layer.Vector('Countries',{
        protocol:  vector_protocol,
        strategies: vector_strategies
    });

    map.addLayer(vector_layer);
    if(!map.getCenter()){
    map.zoomToMaxExtent();

I keep getting thrown this error: mapProjection is null

Does anybody have any ideas what may be going astray?


texts:
I want to import GeoJSON avaliable in the following dump into PostGIS but I am not able to import it. I already imported geousa data available in this public dump into MongoDB smoothly. Let me know is there a way to import this data directly or through MongoDB?

As this data is huge please check following data:

{"geometry": {"type": "Point", "coordinates": [19.056792, 47.490894]}, "type": "Feature",
"id": "SG_1iYphlxn9BSHyGrpv1aXKc_47.490894_19.056792@1308163237",
"properties": {"website": "http://mito.hu", "city": "Budapest",
"name": "Mito Europe", "tags": ["online", "communication", "design",
"branding", "development", "mito"], "country": "HU", "classifiers":
[{"category": "Professional", "type": "Services", "subcategory": "Advertising"}], "href":
"http://api.simplegeo.com/1.0/features/SG_1iYphlxn9BSHyGrpv1aXKc_47.490894_19.056792@1308163237.json",
"address": "N\u00e1dor u. 23.", "owner": "simplegeo", "postcode": "1051"}}

http://s3.amazonaws.com/simplegeo-public/places_dump_20110628.zip


texts: Clipping raster with another raster where selected cell value defined using ArcPy?

texts:
Is there where I can clip a raster with another raster by selecting specific cell value defined?

For example, I want to perform the clip operation using raster "A" as input and using raster "X1" as a clip raster where only interested on cell value = 1 in raster "X1". I have time series of rasters i.e. X1, X2, X3 and want to save the clip results with prefix from clip raster.

Is there way to perform this operation using ArcGIS 10.2 and python?


texts: Adding multiple styles to single layer with pyQGIS?

texts:
On QGIS, I would like to programmatically assign multiple .qml styles to a single layer with python.

I can add styles manually from the Properties dialog on the layer by  selecting Styles->Add... but I cannot add them from a .qml

So the question is two-fold:

How to add programmatically a second, third, etc style for the layer.
If the new styles can be added from  saved .qml's

So far my knowledge is limmited to loading one individual .qml with

layer.loadNamedStyle(path_to_qml_1)


texts: Creating rotated fishnet from grid of evenly spaced points in ArcGIS Desktop?

texts:
I am working with climate models, so I have a point shapefile representing the grid of data points of my climate models. My points are sort of evenly spaced along both directions, the distance being 1,500 m.

This grid is rotated but the rotation angle was not easy to determine, so I skipped this passage. I need to create a fishnet that follows the same rotation as the grid, so that there is one and only one grid point in each square of the fishnet.

How can this be done in ArcGIS?

If I assign the coordinates of my four corners, then they are not included in the fishnet but lay on the border. I want the fishnet to "begin" with an offset of 750 m in each direction with respect to the four corners.

My version is 10.3.1.

This is an image of the grid and a made up fishnet to convey the idea.


texts:
i recently updated QGIS to 2.14 (Essen).
Unfortunately, the very useful plugin << Raster terrain Analysis >> seems to disappear after the update. It was easily accessible by menu Raster, and now is no more. I tried to reinstall it by the plugin manager - it's not even in the list anymore.
The script seems to be correctly in the plugin folder  ( /Applications/QGIS.app/Contents/PlugIns/qgis )

It is deprecated ? How can I reinstall it ?
(i'm on Mac and tried to reinstall QGIS and GDAL framework - nothing).


texts: Mapfish Print rendering incorrect Scalebar even with Geodetic setting

texts:
I'm trying to add a scalebar to my Mapfish Print pdf output but it's always rendering too small.
I have seen a few solutions that suggest using the "geodetic": true setting in my requests, but this makes no difference to the results.
Here's a section of the pdf output demonstrating the issue:

The blue grid marked 75 - 76 is 1 km in size, yet you can see that the scalebar rendered suggests this would be in excess of 1 km if used as the source of truth.
Any ideas how I may be able to get the scalebar to render correctly?
Here's an example request I'm sending to Mapfish Print:
{
    "layout": "A4 landscape",
    "outputFormat": "pdf",
    "geodetic": true,
    "attributes": {
        "map": {
            "projection": "EPSG:3857",
            "dpi": 300,
            "rotation": 0,
            "areaOfInterest": {
                "display": "NONE",
                "area": {
                    "type": "MultiPoint",
                    "coordinates": [
                        [
                            19092823.11,
                            -5325062.01
                        ],
                        [
                            19087166.77,
                            -5329055.85
                        ]
                    ]
                }
            },
            "layers": [
                {
                    "baseURL": "http://tiles-a.data-cdn.linz.govt.nz/services;key=SECRET/tiles/v4/layer=50767/EPSG:3857/{z}/{x}/{y}.png",
                    "type": "OSM"
                }
            ]
        }
    }
}

And the following is my config.yaml:
templates:
  #===========================================================================
  A4 landscape: !template
  #===========================================================================
    reportTemplate: A4-landscape.jrxml
    attributes:
      map: !map
        maxDpi: 300
        width: 802
        height: 520
      overlayLayers: !staticLayers
        default:
          layers:
            - type: "grid"
              numberOfLines: [10, 10]
              labelColor: rgba(0,0,0,0)
              haloColor: rgba(0,0,0,0)
      scalebar: !scalebar
        width: 230
        height: 35
      northArrow: !northArrow
        size: 35
        default:
          graphic: 'north.svg'
    processors:
    - !addBackgroundLayers
      inputMapper:
          overlayLayers: staticLayers
          map: map
    - !reportBuilder # compile all reports in current directory
      directory: '.'
    - !createMap {}
    - !createScalebar {}
    - !createNorthArrow {}


texts: Converting ASCII file to DEM using SciPy?

texts:
ArcGIS does not employ a pure linear interpolation algorithm, so I turned to python and specifically scipy to perform the interpolation.

I found out I can create an ASCII file containing the interpolated values and pass that to ArcMap in order to come up with a DEM. I successfully performed the interpolation, then created the ASCII per the instruction provided in the ESRI ASCII Raster format documentation and converted the ASCII to a raster using ArcMap's corresponding tool.

I first tried importing random values in order to check that this procedure actually works and it worked as expected. But when I pass the interpolated values, although the conversion tool does not raise any errors, the DEM makes no sense.

Below is the python code I used and an image of the DEM.

Can anyone solve any issues or recommend a workaround or another way to do it?

import sys
import numpy as np
from scipy.interpolate import griddata

np.set_printoptions(formatter={'float_kind':'{:f}'.format})
x,y,z = np.loadtxt(sys.argv[1],delimiter=',',usecols=(1,2,3),unpack=True)
pixel = np.float(sys.argv[3])
xmin = min(x)
ymin = min(y)
xmax = max(x)
ymax = max(y)
nrows = (ymax - ymin)/pixel
ncols = (xmax - xmin)/pixel
dir_txt = sys.argv[1].replace('.txt', '')

xi,yi = np.mgrid[xmin:xmax:complex(0, ncols), ymin:ymax:complex(0, nrows)]
zi = griddata((x, y), z, (xi, yi), fill_value=-9999, method=sys.argv[2])

# ASCII file header
header = "ncols %d \nnrows %d\nxllcorner %f\nyllcorner %f\ncellsize %f\nnodata_value -9999" % (int(ncols), int(nrows), float(min(x)), float(min(y)), pixel)
np.savetxt("out.txt", zi, fmt='%f', header=header, comments='', delimiter=' ', newline='\n')


texts: Monster Spatial database server suggestions

texts:
I'm going to setup a open Mapserver w/ Openlayers, Geojango, Postgis using Postgresql. And I need some suggestions from people that have been shopping for what I need longer than I have.

What would be the best configuration for such a db server and from where?

Im a big mac fan and was thinking of a Xserve w/ 42g ram and two quad cores (2.66). What does postgresql like best as far as server software/hardware? Can it fully utilize multi cores?

Im on somewhat of a tight budget so price is a big factor. But what would be the dream machine and its cheaper inequivalent?

sockets; i'd probably have 20 to 30 ppl at the lease almost 100 at most, maybe more. Serving up more than 100 layers with some common canned geoprocessing type stuff going on.


texts: Accessing GeoServer WFS services in ArcGIS Online map viewer?

texts:
Has anyone was able to add GeoServer WFS service into ArcGIS Online map viewer? By this I mean having an ArcGIS Online account (either free or for Organizations), logging in and adding layers from the web.

I have no problems adding a WMS service such as this one
http://demo.opengeo.org/geoserver/wfs?service=wfs&version=1.0.0&request=getfeature&typename=topp%3astates. It is being retrieved and drawn nicely in ArcGIS Online map viewer.

What about the WFS services? When trying to add this service (which should be a WFS?) http://suite.opengeo.org/geoserver/wfs?service=wfs&version=1.1.0&request=GetFeature&typename=usa:states&featureid=states.23, I am getting only images (somehow "redirected" to WMS maybe?). OGC things are not really my strongest skills.

Is it possible to get raw vector data (served by GeoServer) instead of map images into the ArcGIS Online map viewer? Any workarounds to get it from GeoServer via some other source into ArcGIS Online?


texts: I am having a problem with related table in collector for ArcGIS. I have a feature class with a related table defined by a relationship class of type simple, 1:1 cardinality, and forward notification. I have also activated the pop-up for both the feature class and the table in ArcGIS online. When I enter data in Collector in feature class, the related table appears and I can enter data into it too. After pressing "Done", when I try to access the related table data tapping "view", I get an error saying that the data could not be found. When I access the related table in ArcGIS online, I can see that the data I entered are stored in the table but as soon as I try to access it with pop-up, I get the same error. Can anyone tell me, what I am missing here or what I am doing wrong?


texts:
I understand how to capitalize first word using !CITY_NAME!.capitalize().  However, if I have a string in a record and I want to capitalize the fifth word, what is the code for that?

Ex: 3 Boundary in fifth column How do I capitalize column?

I can use: 'column'.replace('c','C') but is there a way to call the capitalize function and not replace to do it?


texts: How to display georeferenced image in polygon in OpenLayers 3?

texts:
I have a GeoJSON polygon whose coordinates represent the four corners of an image I want to display on an OpenLayers map.

{
      "coordinates": [
          [
              [16.15799, -75.81512],
              [22.09792, -76.89212],
              [17.24105, -78.28423],
              [11.10162, -77.09659],
              [16.15799, -75.81512]
          ]
      ],
      "type": "Polygon"
}

I see plenty of examples using ol.source.ImageStatic to overlay an image, using a [minx, miny, maxx, maxy] extent but this results in an image displayed aligned to the top of the canvas.

  const extent = [ coordinates[0][0][0],
                   coordinates[0][0][1],
                   coordinates[0][2][0],
                   coordinates[0][2][1]];

  const projection = map.getView().getProjection();
  const imageExtent = ol.proj.transformExtent(extent, 'EPSG:4326', projection);
  const layer = new ol.layer.Image({
    source: new ol.source.ImageStatic({
      url: image.properties.preview_image,
      imageExtent,
      projection,
    }),
  });

What approach should i take to fit the image within the polygon?

I originally tried taking the bounding extents from the polygon, as suggested in Ivan's answer below, but this results in the image being overlaid on the polygon at the wrong orientation, between two points that aren't on its edges:

const format = new ol.format.GeoJSON();
const feature = format.readFeature(image);
const extent = feature.getGeometry().getExtent();

I've read threads such as this one which suggests manipulating the canvas with the precompose and postcompose events of the map layer, but this seems pretty deep into the API just to display a georeferenced image:
https://groups.google.com/forum/#!msg/ol3-dev/bOrpFyRIixc/I1BZMOmwBQAJ


texts: How does the IdentifyTask handle proxy/cross-origin issues?

texts:
Background:

I'm trying to perform an Identify against an ArcGIS Server layer, without using the ArcGIS Server JS API, and am running into issues which don't occur when using the ArcGIS Server JS API.

Download the Display identify results in popup sample to your local web server
Click on a parcel, and note that the Identify task runs against the sample server using these parameters. The results are correctly returned to the browser.
In a simple HTML page with jQuery enabled, I'm trying to run the equivalent query using a direct AJAX call to the layer's Identify task on the ArcGIS Server REST API:

  var geometry = decodeURIComponent('%7B%22x%22%3A' + -9270290 + '%2C%22y%22%3A' + 5247184 + '%7D');

  var data = {
    "sr": 102100,
    "tolerance": 3,
    "maxAllowableOffset": 0.1,
    "returnGeometry": true,
    "imageDisplay": "1920,938,96",
    "mapExtent": "9270703,5246922,-9270168,5247482",
    "geometry": geometry,
    "geometryType": "esriGeometryPoint",
    "f": "json",
    "layers": "all:0,2"
  }

  $.ajax({
    dataType: "json",
    type: 'GET',
    url: "https://sampleserver3.arcgisonline.com/ArcGIS/rest/services/BloomfieldHillsMichigan/Parcels/MapServer/identify",
    data: data,
    success: function(result){
      console.log("success");
    },
    error: function(error){
      console.log("error");
    }
  });

The resulting URL works in the browser, but the AJAX call fires the error function, with a cross-origin error message: No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://localhost' is therefore not allowed access:

My question:

How is the IdentifyTask.execute function (from in the ArcGIS Server JS API) able to avoid this cross-origin error? And/or how can I emulate that behaviour using a direct AJAX query?


texts: How to create raster from table containing cell values in PostGIS?

texts:
I am having a table with integer values. These values are in between 0 and 2. I want to create a raster out of that table with integers as cell values. Or else how to create polygons out of that table with each polygon consisting of same integer value.

That table is not geometry. So, is it possible to check whether the particular pixel value is present in the raster or not. If I can get the this I can use ST_SetValues for assigning new pixel value.


texts: Remove multiple fields from multiple shapefiles using a python script

texts:
I'm using the python console of QGIS 2.8.1-Wien.
I have a shapefile with many fields in the table that I want to remove them.
I tried to use this geoalgoritm as:

processing.runalg('qgis:deletecolumn', myshapefile.shp, "myfield", myoutput.shp)

This works perfectly and the field I want is removed.
My question is: if I want to remove more than one fields how should I write this formula? Or there are any other way to write a python script?

I used your 1st method, just for one shapefile as:

output_0=processing.runalg('qgis:deletecolumn', "D:/PATRICIA/ECMWF/Test/TPrecip_20140628_VERAO_join9.shp", "Lat_dec", None)
output_1=processing.runalg('qgis:deletecolumn', output_0, "Lon_dec", None)
output_2=processing.runalg('qgis:deletecolumn', output_1, "Altitude", "D:/PATRICIA/ECMWF/Test/TPrecip_20140628_VERAO_join12.shp")

The result was:"Error: Wrong parameter value: None"

I tried to use the 2nd method for 3 similar shapefiles (with same fields) activated in the table of contentes but just occured de removal of the fields in the 1st shapefile :

layer = qgis.utils.iface.activeLayer()
res = layer.dataProvider().deleteAttributes([1, 2, 3, 4, 5, 6, 7, 8])
layer.updateFields()

In this case it's impossible to use 8000 shapefiles activated in map canvas.

Any help?


texts: Why are layers not overlapping to each other inspite of having same CRS?

texts: I am new to QGIS and working with QGIS 1.7.3 version, I have my all layers in WGS 84 CRS but even though they are not overlapping to each other, What should I do ?


texts: Opening multiple text files into a single .tab file in MapInfo

texts:
I am trying to open multiple xyz text files into mapinfo to create a table where i can create a raster image of multibeam bathymetry.
Is there a script which will open each .txt file and merge them all into a single .tab file. And even create a raster image from this?


texts: I have four fields in my shapefile. I would like to fill only one of them (cell of attribute table) and automatically after to write one cell the others one will be filled by conditional clause. I need that this conditional has to keep there, filling the value attribute, without necessary to run filed calculator (automatically).


texts: Copy feature class from multiple gdb to a single gdb

texts:
I need to write a Python script that is copying. I want to copy a common feature class of these gdb to a single gdb

import os
import arcpy

arcpy.env.workspace = ("D:/Planchas")

workspaces = arcpy.ListWorkspaces("", "Access")

print(workspaces)

for workspace in workspaces:
    featureSet = arcpy.ListDatasets("Toponimos")
    for fs in featureSet:
        featureClasses = arcpy.ListFeatureClasses("Orografia", "Point", fs)
        for fc in featureClasses:
            arcpy.CopyFeatures_management(
                    fc,
                    os.path.join("C:/Users/Prueba",
                    os.path.splitext(fc)[0]))

When I try to run the code there is no error but the process is not carried out


texts: How to prevent QGIS (MMQGIS plugin) from interpreting numbers as strings

texts:
I'm using the MMQGis plugin to "Join Attributes From CSV" into a shapefile and am having trouble getting QGIS to interpret some numbers as numbers (rather than strings which is what its doing).

This is becoming a real problem because the values do not behave like they should when trying to order them, select above a range, etc.

I've seen two somewhat similar questions to this but the answer to this question instructs the user to import the CSV as a layer. But the CSV data I'm trying to import is not geocoded and I'm using values (things like zip code, etc) to join in the CSVs, so I don't think that question's answers might help here.


texts:
I have been using geoserver from quite sometime now and publishing layers and consuming services. Recently I have started facing the issue mentioned below when I am trying to preview the layer in geoserver using openlayers, but renders the layer using geojson. I am unable to resolve this. Kindly provide pos

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://x.x.x.x:8080/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException code="ForbiddenFormat">
      Creating maps using application/openlayers is not allowed
</ServiceException></ServiceExceptionReport>


texts:
I am using ArcMap 10.8.
I add the UAV data for raster data(.tif).
But in my display, the data only display a particular scale(1:600).
I checked my scale range. But also set "show layer at all scales".
What can I do in this situation?


texts: Drill trace line (arrow) from a point (collar) along a given azimuth with a given length

texts: I'm trying to create a diamond drill surface trace line (arrow) from a point (collar) along a given azimuth with a given length from a csv file. This actually stems from a problem of projecting a surface trace from a drill string in Diamond Drilling. I've used trig to define the surface trace on a flat surface but would love to tie this to a DEM in the future!


texts: Creating multiple donuts with ascending hole radius (Multiple Ring Buffer)

texts:
I want to use Multiple Ring Buffer (MRB) tool to create multiple donuts, but when I use NONE dissolve and outside polygon it cuts only input buffer from all created and I want it to cut first ring from the second one created.
For better understanding:
I create 3 buffers - 1km, 5km, 10km
First is continuous from starting point and I use it as input for MRB tool.
2nd and 3rd created via MRB:
2nd is what I want because tool dissolves and is outside it.
3rd - I want to have hole with 5km radius=cuted 2nd buffer not 1st.

I know I can do it creating one after another but maybe there is way to make it in one operation?
I assume that there is small (?) change of python code needed but I don't know where and how :(


texts: Iterative position within Field Calculator/Python

texts:
I'm not great with Python or programming in general, so I'm not even sure how to ask this question.

Basically what I've got is in the left column, and what I want is in the right column.

UniqueID    POSITION
   a              1
   a              2
   a              3
   a              4
   a              5
   b              1
   c              1
   c              2
   d              1
   e              1

I guess what I need to do is loop through the sorted UniqueID field and increment the Position field UNTIL the UniqueID changes, in which case the Position should reset to 1.

NFI how to make this actually work.  I've done a lot of searching and have found some promising-looking code, but I just don't have the know-how to make it fit my problem.  Sorry about my brain.  My script right now looks like "Baby's first adventures with arcpy," and since this is slightly more complicated than a simple arcpy call, it's stopped me in my tracks.

I guess I'm looking for something that works as part of a CalculateField operation?


texts: Plotting shapefile in R: how to delete selected borders?

texts:
I'm plotting a shapefile in R and would like to delete certain borders to "group" shapes together. For example, in the plot below, how would I delete the borders between the three blue counties and between the red counties?

library(maptools)

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
 IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))

COLS <- rep(2,100)

COLS[c(9,31,82)] <- 4

plot(xx, axes=TRUE, las=1, col=COLS)


texts: Add corresponding MODIS LAI band to Sentinel-1 log(VH/VV) data in Google Earth Engine

texts:
I am trying to determine the correlation between Sentinel-1 log(VH/VV) data and the corresponding MODIS 4 day LAI product.  For the first step, I am trying to add the corresponding LAI data as a band to the Sentinel-1 log(VH/VV) collection.  My code so far is provided.  I am getting the errors shown below.

var point = ee.Geometry.Point(-96.47664, 41.16506);
var proj = ee.Projection('EPSG:32614');

var reproject = function(image) {
  var reprojected = image.reproject(proj, null, 1000);
  return reprojected};

var LogRatio = function(img) {
  var logVV = img.select('VV');
  var logVH = img.select('VH');
  var logratio = logVH.subtract(logVV).rename('logVH/VV');
  var logratio_float = logratio.toFloat();
  return img.addBands(logratio_float)};

var S1 = ee.ImageCollection('COPERNICUS/S1_GRD');
var S1_TempFilt = S1.filterDate('2017-01-01', '2020-02-20');
var S1_TempFilt_Resamp = S1_TempFilt.map(reproject);
var S1_TempSpatFilt_Resamp = S1_TempFilt_Resamp.filterBounds(point);
var S1_TempSpaceFiltResamp_W_LogCR = S1_TempSpatFilt_Resamp.map(LogRatio);
var S1_Log_VH_VV = S1_TempSpaceFiltResamp_W_LogCR.select('logVH/VV');
print('S1_Log_VH_VV',S1_Log_VH_VV);

var MODIS = ee.ImageCollection('MODIS/006/MCD15A3H');
var MODIS_TempFilt = MODIS.filterDate('2017-01-01', '2020-02-20');
var MODIS_TempFilt_Resamp =  MODIS_TempFilt.map(reproject);
var MODIS_TempSpatFilt_Resamp = MODIS_TempFilt_Resamp.filterBounds(point);
var M_LAI = MODIS_TempSpatFilt_Resamp.select('Lai');
print(M_LAI);

var Find_MODIS_LAI = function(image2) {
  var S1date = ee.Date(image2.get('system:time_start'));
  var LAI_start_date_min = S1date.advance(-1,"day");
  var LAI_start_date_max = S1date.advance(2,"day");
  var MODIS_LAI = M_LAI.filterDate(LAI_start_date_min, LAI_start_date_max);
  return MODIS_LAI};

var MODIS_LAI = S1_Log_VH_VV.map(Find_MODIS_LAI);
print('MODIS_LAI', MODIS_LAI);

var LAI_Band = MODIS_LAI.select('Lai');
print('LAI_Band',LAI_Band);
var LogCR_LAI = S1_Log_VH_VV.addBands(LAI_Band);

Error:

ImageCollection (Error)
Error in map(ID=S1A_IW_GRDH_1SDV_20170108T002112_20170108T002137_014733_017FB8_9BC3):
Image.select, argument 'input': Invalid type. Expected: Image<unknown bands>. Actual: ImageCollection.

S1_Log_VH_VV.addBands is not a function


texts: Cannot add Leaflet markers to map iterating through CSV

texts:
I have data from a CSV that I want to turn into markers using the latest version of LeafletJS. I can successfully turn the CSV into an array of js objects and have confirmed that the latlng information is correct and being properly created. However, the markers do not place on the map.

// create map object
var mymap = L.map('mapid');

// fetch tiles from Mapbox API
L.tileLayer('https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token={accessToken}', {
    attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © Mapbox',
    maxZoom: 18,
    id: 'mapbox/streets-v11',
    tileSize: 512,
    zoomOffset: -1,
    accessToken: [removed for privacy]
}).addTo(mymap);

// set view on page load
$(document).ready(function() {
    // [40.47, -79.98]: LatLng for Allegheny County
    // 12: zoom parameter
    mymap.setView([40.47, -79.98], 12)

    $.get('../data/output.csv',function(data) {
        var results = $.csv.toObjects(data);
        console.log("printing")
        console.log(results);
        results.forEach(result => {
            var latlng = L.latLng(parseFloat(result.Latitude),parseFloat(result.Longitude));
            var marker1 = L.marker(latlng).addTo(mymap);
        });
     });
});

Does anyone know why the markers aren't being created?


texts: Raster Version of Spatial Join (one to many)

texts:
I have two raster layers that are identical in shape (they overlap exactly) and I need to perform the raster equivalent of a Spatial Join (one to many). My files need to stay in raster format otherwise I would just do everything in vector. I'm using ArcMap 10.1.

My "buffrast" raster contains the land cover value for each pixel within a circular buffer of 150m around a bunch of pointa. My "rr1" raster contains the identities for individual 150m radius circles around various points (e.g. NR0602), but has no land cover value.

I need to have a final raster whose attribute table lists the land cover value for each pixel within a given 150m radius circle. That is, I need to have the circle's identity (e.g. NR0602) listed for each pixel along with its land cover value. I've already tried raster mosaic, raster calculator (combine) but I think my syntax must be incorrect.

I am completely new to raster analysis so please be as explicit as possible in any advice you might have.

Thanks!
alyssa


texts: Geoserver Point-to-Line Aggregation possible?

texts:
I'm trying to use Geoserver to serve a WMS displaying lines instead of points, but without changing the underlying data structure. Currently I'm displaying the data as simple points, which works fine.

Is it possible to use geometry / rendering transformations to connect the points to lines based on the timestamp?
Or is there a plugin that would allow me to aggregate the points to lines? For example by creating a lineString from the point's coordinates.

Environment:

Geoserver 2.11.0
ElasticSearch 5.3.2
ElasticGeo 2.11.1-RC1 (Plugin    that allows Geoserver to use
ElasticSearch as Datastore)

Elasticsearch data type:

lat / lon mapped as GEO_POINT
timestamp mapped as Date string or UNIX-based


texts: Uncompressed size of ASTER DEM tile is 0 Bytes

texts:
I have downloaded ASTER DEM tiles(n=400), I have to merge/mosaic them. While doing so I got a problem with one tile which does not merge but displays correctly.

I have downloaded it again but still I get the same problem. How can I mosaic it with others. I have seen such a question on this forum but it has not been answered. Link. The only properties it has different with others is that its uncompressed size is 0 bytes, no columns nor resolution. I have realized that the bad tile, one is able to view the attribute table while the rest one can not. The tile number is ASTGTM2_N30E041_dem.tif


texts:
Trying to make a couple databases which contain everything that needs to be mapped for a few projects. For best conformity it would be great to keep style embedded within the database to enable a simple save as but the only way I can think of it to .qml all the different styles. As there will be >20 layers all in, don't really want a different .qml for each layer, takes up space and less easy!

But the symbology export doesn't ever seem to keep the symbology.

Can it be done?

Can save as different file type but wouldn't know which one is best.


texts: I have a line layer (rivers) and a polygon layer (forests). Both files have the same SRS. I want in a new layer clear all the rivers that do not flow through the forest. When I try clip tool, I have this error


texts:
I have been trying to load a shapefile available at https://www.elections.on.ca/en/voting-in-ontario/electoral-district-shapefiles/open-use-data-product-licence-agreement/download-shapefiles.html to SQL Server 2014 in the EPSG: 4326 formats.

I tried this using the Shape2SQL tool available at https://www.sharpgis.net/page/SQL-Server-2008-Spatial-Tools

I checked the "Planar Geometry" option and left the "Set SRID" field unchecked. If I check the "Set SRID" field and enter 4326 in the field, nothing gets uploaded to SQL Server. Hence, I left it unchecked.

The loaded coordinates look like this and are NOT in LAT LONG form:

ID  ED_ID ENGLISH_NA       FRENCH_NAM       SHAPE_area    SHAPE_len  geom
1   8     Brampton Centre  Brampton-Centre  4.444118E+07  32761.97   0x000000000104B6010000A...

After uploading to the database, when I checked the SRID for the imported table in SQL Server using the command:

select distinct SP_GEOMETRY.STSrid from dbo.MYTABLE

It returns "0".

Please help with following questions:

Is the shapefile getting loaded into the SQLServer correctly?
How do I convert the loaded data to EPSG: 4326?


texts:
I would like to observe things like maps.on('moveend') to update zoom level or center values. Currently I have a "construct" made up of computed params and watcher functions like this:
export default {
  [...],
  // Computed values are saved / read from vuex store
  computed: {
    currentLatLon: {
      get () {
        const obj = {
          lat: this.$store.getters.getCurrentLat,
          lon: this.$store.getters.getCurrentLon
        }
        return obj
      }
    },
    moveEnd: {
      get () {
        const val = this.olMap !== null
          ? this.olMap.on('moveend', () => {
              console.log('moveEnd; do something... getCenter()...'))
            }
          : false
        return val
      }
    }
  },
  watch: {
    currentLatLon (value) {
      console.log('watcher currentLatLon invoked: value =', value)
      // this.zoomToCurrentLocation()
    },
    moveEnd (value) {
      console.log('moveEnd: value =', value)
    }
  }
}

codesandbox with minimal code:
https://codesandbox.io/s/openlayers-getzoom-74gnw?file=/src/components/MapContainer.vue
It works so far. But I am unsure, because it seems to be "double the trouble".
Is there a better way to observe Openlayers events in Vue?


texts:
I have created a regular grid of (square) polygons using Processing->Create grid.
I would now like to iterate over my grid squares and clip an underlying road network within each grid square. I can do this from the GUI but was wondering how to do it using pyQgis (or console).

With l1 the layer of lines and l2 the grid
My attempts have been along the lines of

iter=l2.getFeatures()
for f in iter:
    processing.runalg("qgis:clip", l1, f.attributes(), None)

But this fails with a wrong parameter... The third argument should be a polygon or vector. How do I convert?

Where am I going wrong?


texts: Subscript out of bounds error while mosaicking TIFFs together

texts:
I have just downloaded NDVI imagery from Google Earth Engine. The image was too big so it was divided into about 35 .tif files. Each tile has approximately 950 bands of NDVI values.
I am now trying to stitch the images back together in R. I have tried two different methods:
The first method which uses the raster package is really slow:
library(raster)

f <- list.files(path = "/scratch/MODIS_NDVI_TIFS", pattern = ".tif$", full.names = TRUE)

rl <- lapply(f, raster::brick)

m <- do.call(raster::mosaic, rl)

writeRaster(m,"APL_NDVI_Combined.tif")

The second method I found discussed here which I modified for my use as such:
library(gdalUtils)
library(rgdal)
library(raster)

f <- list.files(path = "/scratch/MODIS_NDVI_TIFS", pattern = ".tif$", full.names = TRUE)

e <- extent(113.338953078, 153.569469029,-43.6345972634, -10.6681857235)

template <- raster(e)
proj4string(template) <- '+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs '

writeRaster(template, file="Combine_NDVI.tif", format="GTiff",overwrite=TRUE)

mosaic_rasters(gdalfile=f,dst_dataset="Combine_NDVI.tif",of="GTiff")

However, I get an error in during the moasic_rasters call which states:
Error in strsplit(gsub(grep(cmd_output, pattern = "Size is ", value = TRUE),  :
subscript out of bounds
Calls: mosaic_rasters -> gdalinfo
Execution halted

Any idea what is going on here? Are there better suggestions for combining these rasters together?


texts:
i have an hdf file with x,y (pixel/line) datasets as viewed in HDFView. When i convert it to tif file this dataset does not appear. How should i make gdal read the georeferencing dataset in order for me to project it properly?
these are the subdatasets:

these are the datasets as shown by HDFView.


texts: Bring a list of unsorted linestrings in order so that they form a path from origin to destination

texts:
I have a table of linestrings together with some attributes that, when displayed all at the same time, represent a connected path from an origin to a destination (each linestring has one predecessor, except for the first one, and one successor, except for the last one). However, the feature order in my table is randomized.

For another application, I need the features in the table to be ordered from origin to destination (or vice versa), so that for each row/linestring in the table, the preceding row represents the current row's predecessor linestring and the succeeding row represents the current row's successor linestring.

I am in control over the data, so writing them out in any format suitable for this operation should be possible (csv, geojson, shapefile, ...).

I have QGIS, ArcMap, PostGIS and Python easily available but any other solution will also be fine.

As an example that everybody can easily visualize, take the following GeoJSON:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "someAttribute": 2 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88175, 52.294 ], [ 6.88198, 52.29404 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 8 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88945, 52.29465 ], [ 6.88955, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 7 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88934, 52.29465 ], [ 6.88945, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 4 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88294, 52.29418 ], [ 6.8857, 52.29448 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 3 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88198, 52.29404 ], [ 6.88294, 52.29418 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 5 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.8857, 52.29448 ], [ 6.88836, 52.29462 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 9 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88955, 52.29465 ], [ 6.88963, 52.29466 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 6 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88836, 52.29462 ], [ 6.88934, 52.29465 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 11 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89216, 52.29463 ], [ 6.89285, 52.2946 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 12 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89285, 52.2946 ], [ 6.89528, 52.29444 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 1 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88108, 52.29389 ], [ 6.88175, 52.294 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 13 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89528, 52.29444 ], [ 6.89794, 52.29414 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 10 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.88963, 52.29466 ], [ 6.89216, 52.29463 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 14 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.89794, 52.29414 ], [ 6.90067, 52.2937 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 16 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.9034, 52.29311 ], [ 6.90369, 52.29304 ] ] } },
{ "type": "Feature", "properties": { "someAttribute": 15 }, "geometry": { "type": "LineString", "coordinates": [ [ 6.90067, 52.2937 ], [ 6.9034, 52.29311 ] ] } }
]
}

For clarification, if you visualize this (e.g. with QGIS), you get the following picture:

The single linestrings form a connected path from 1 to 16, even though they are unordered in the data source. What I want is to sort my data source accordingly. The enumeration in someAttribute of course is only for clarification and not given in the real data, so simply sorting by attribute is not an option ;-)


texts: Is it possible to set symbol layers always on top?

texts:
Is it possible to set symbols always on top, also if the layer with the symbol is drawn under another layer?

Example:
My layer order is like that:

Streets (top)
Buildings
Landuse (bottom)

The Buildings layer have centroid-symbols. But these are drawn underneath the streets layers: