texts: Google earth engine: How to get yearly DN value by satellites from DMSP-OLS nighttime satellite data

texts:
I would like to get the panel data of the nighttime light in Japan and check the correlation between other economic statistics. In order to do that, I have to get yearly digital number value from each satellite, F10, F12, F14, F15, F16, and F18, and intercalibrate it as Elvidge et al(2009) did.
Here is the code in which I could get the DN value in 1994, the year which F10 and F12 exist at the same time. but I can't identify which satellite are the data from.
So, the question is below:
How can I check which satellite are the data from?
If it is impossible to make sure of it, How to get DN value from a certain satellite?
var dataset = ee.ImageCollection('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS').select('stable_lights');
var v1994 = dataset.filter(ee.Filter.date('1994-01-01', '1994-12-31')).median();
var nighttimeLightsVis = {min: 0.0, max: 60.0};
Map.addLayer(v1994.clip(jpn), nighttimeLightsVis, "Nighttime_1994");
var reducer1 = ee.Reducer.mean();
var reducers = reducer1.combine({reducer2: ee.Reducer.median(), sharedInputs: true})
function ntl(image) {var pas = image.reduceRegions({  collection: jpn,  reducer: reducers,  scale: 500,});return pas}
function export_table (table, description) {var em = Export.table.toDrive({ collection: ee.FeatureCollection(table), description: description, folder:'GIS', selectors:(["CITY_ENG","JCODE","mean"]),});return em}
var nl_v1994 = export_table(ntl(v1994), 'NL_v1994Ejpn');

*I got the border data, "jpn" in this code, from this website
*I tried to use the GeoTIFF image which downloaded from this website, but I got an error: Image.reduceRegions: Unable to find a crs.


texts: Landsat 5 imagery not lining up with Landsat 8 imagery

texts:
My Landsat 5 and Landsat 8 imagery are not lining up in ArcMap (also tried QGIS). (See attached images) L8 seems to be ok and lines up with the shapefile where as with L5, the shapefile is farther inland.

The data frame as well as both sets of imagery are in WGS_1984_UTM_Zone_21N.

The shapefile has been projected to the same.

*edit: Sorry for the lack of details :)
The 2017 image (and all L8 images I downloaded from EarthExplorer) line up correctly with the polygon (mangrove restoration site). All L5 images I'm using (2007-2011) place the polygon too far inland.

The polygon was projected to the same crs as the landsat imagery (WGS_1984_UTM_Zone_21N).

Any ideas why the L5 imagery is off? Could it be with the imagery itself? Now that I think about it, I uploaded the shapefile in EarthExplorer, and I do remember one scene that I looked at and the polygon was too far inland.


texts: How to show other_tags columns separately?

texts: I have many columns in an attribute of a file downloaded from OSM. But QGIS only shows three columns including one named "other_tags" that includes all my other columns and their value. How do I show all those columns separately?


texts: How to split polygons with another polygon and inheritage data?

texts:
Let's say that we have two tables with polygons representing areas. Table A and Table B. Geometries from table A represents both the valid and invalid areas. And the geometries in table B  represents the invalid ones.

For example, lets say that the geometries from table A represent, harvesting fields, and table B represents roads, shacks and generally areas that the harvesting is impossible for w/e reason (as polygons). Table B represents some kind of mask that's has being digitized differently.

The problem arises that the mask, can pass through many fields. So I want to clip that mask (Table B), with the field polygon (Table A) to create a new set (Table C) consisting with the  geometries that represents the masks 'per field' along with the name of the field it was cut from.

This can be done by

create tableC as (select st_intersection(a.geom,b.geom) geom from tableA as a, tableB as b);

But I have troubles on how to select the tableA.name as well and transferring it to tableC.

the selection part of an ideal querry would be like this:

create tableC as (select st_intersection(a.geom,g.geom) geom, tableA.name as name from .... );

A sample defintion of table C is the following:

Column     |   Type      |
gid        | int         |
geom       | geometry    |  <- st_intersetion(a.geom,b.geom)
name       |varchar(255) |  <- name of the feature of table A which took part in the st_intersection

Table A is defined as  following:

Column     |   Type      |
gid        | int         |
name       | varchar(255)|
geom       | geometry    |

And Table B as follows:

Column     |   Type      |
gid        | int         |
geom       | geometry    |

And for a visual representation :


texts: Is it possible to access javascript variables that are in an OBJECT tag?

texts:
I want to inject a whole html page into a "main" page. I succeeded in doing this, but I still need to access variables of the injected page, from within the main page. Here is an example:

I have my index page, where I injected a GIS page.

I want to access the map variable (created in "gis.html") from within my index page, in order to control my GIS app from the main (index) page.

index.html

        <title>Index</title>
        <meta charset="UTF-8">
        <meta name="viewport" content="width=device-width">

            <button id="zoomOut" onclick="zoomOut()">ZOOM OUT</button>
            <button id="zoomIn" onclick="zoomIn()">ZOOM IN</button>

gis.html

        <title>GIS</title>
        <meta charset="UTF-8">
        <meta name="viewport" content="width=device-width">
        <link rel="shortcut icon" href="#">

        <link href="https://cdnjs.cloudflare.com/ajax/libs/ol3/3.16.0/ol.css" rel="stylesheet" type="text/css">

Here are the pages that I am talking about:

index
gis

EDIT 1: DIFFERENT DOMAINS ---> IMPORTANT

The two html pages are supposed to be on different domains. For example, "gis.html" is the GIS Web App that I developed. The "index.html" would be a webpage of a client that wants to integrate a small web app on his site.


texts: Partial loss of projection data converting to MapInfo File with OGR

texts:
I'm trying to convert some data from ESRI Shapefile to MapInfo TAB format using ogr2ogr. The data is projected using the British National Grid (EPSG:27700). The output file is valid, and the geometry and attribute data are OK, but the projection is incorrect. This is best illustrated with an example.

Start with a GeoJSON file containing a single point:

{
    "type": "FeatureCollection",
    "features": [{
        "type": "Feature",
        "properties": { "FID": 0 },
        "geometry":
        {
            "type": "Point",
            "coordinates": [ 393800.0, 212700.0 ]
        }
    }]
}

If I convert to an ESRI Shapefile (assigning a projection, as GeoJSON has none), everything works:

$ ogr2ogr test.shp test.geojson -a_srs EPSG:27700
$ gdalsrsinfo -o proj4 test.shp
'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs '

QGIS expands +datum=OSGB36 to a 7 parameter transformation from WGS84, as so:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

Next, I try to convert to MapInfo:

$ ogr2ogr test2.tab test.geojson -a_srs EPSG:27700 -f 'MapInfo File'
$ gdalsrsinfo -o proj4 test2.tab
'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +units=m +no_defs '

The result is a 3 parameter transformation. The actual geometry is unmodified, and I can manually redefine the projection in QGIS, but this is not good.

Next, I try to reproject to a 3 parameter transform...

$ ogr2ogr test3.tab test.geojson -s_srs EPSG:27700 -t_srs '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +units=m +no_defs ' -f 'MapInfo File'
$ gdalsrsinfo -o proj4 test3.tab
'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

It changed again! Finally, if I try to reproject once more:

$ ogr2ogr test4.tab test.geojson -s_srs EPSG:27700 -t_srs '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ' -f 'MapInfo File'
$ gdalsrsinfo -o proj4 test4.tab
'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

This works - the geometry has been modified, but at least it is what it says it is.

But that's not what I want! How to I get the MapInfo file to store the 7 parameter transform?

This issue seems to be related to these questions, but as far as I can see they don't provide an answer. From what I can tell, either OGR is doing something wrong, or the MapInfo format simply can't handle this information.
QGIS and MapInfo Projections - EPSG 27700 (British National Grid)
Question about BNG coordinate systems


texts: Methods behind qgis:simplifygeometries and qgis:smoothgeometry?

texts:
I would like to know on which method (e.g. Douglas-Peucker line simplification) the modules qgis:simplifygeometries and qgis:smoothgeometry are built on.

I cannot find further information on these modules.


texts:
I am creating a Python/PyQGIS tool in QGIS 2.18 which combines a few functions:

Extract several attributes from source layer
Redefine a coordinate system
Optional: clip the source layer with a clipping layer

The goal is to run this tool as a batch processing over multiple layers.

The tool is working now, but results sometimes in empty layer because the source layer does not have features within the clipping layer. In that case I want to skip the output (no output). Does anyone know how to built that in Python?

This is my current code:

##BGT_to_shape_nieuw2=name
##input=vector
##clip=optional vector
##output=output vector

from qgis.core import *
from qgis.core import QgsFeature

layer = processing.getObject(input)
select = QgsVectorLayer.selectByExpression(layer, 'eindRegistratie is NULL')
num = layer.selectedFeatureCount()
if num > 0:
    if clip is not None:
        outputs_QGISEXTRACTBYEXPRESSION_1=processing.runalg('qgis:extractbyexpression', input,'eindRegistratie is NULL',None)
        outputs_QGISDEFINECURRENTPROJECTION_1=processing.runalg('qgis:definecurrentprojection', outputs_QGISEXTRACTBYEXPRESSION_1['OUTPUT'],'EPSG:28992')
        outputs_QGISCLIP_1=processing.runalg('qgis:clip', outputs_QGISDEFINECURRENTPROJECTION_1['OUTPUT'],clip,None)

    else:
        outputs_QGISEXTRACTBYEXPRESSION_1=processing.runalg('qgis:extractbyexpression', input,'eindRegistratie is NULL',output)
        outputs_QGISDEFINECURRENTPROJECTION_1=processing.runalg('qgis:definecurrentprojection', outputs_QGISEXTRACTBYEXPRESSION_1['OUTPUT'],'EPSG:28992')


texts: Setting values for palette in Google Earth Engine

texts:
I am very new to GEE. I have been trying to figure out how to assign different colours to different NDVI values. Let's say, in following code, I want to assign 'red' for NDVI value of 0.2-0.4, 'orange' for 0.4-0.6, 'yellow' for 0.6-0.8, and 'green' for 0.8-1. When I use following code, I feel that I do not have full control over the colours that are being assigned to the NDVI values. The colours I get are gradual based on the NDVI value at certain pixel. What I am trying to do is displaying all pixels that fall in the respective NDVI range in the same colour (something that is done in symbology of ArcMap or QGIS).
link to the code: https://code.earthengine.google.com/1f92036224ea8e7e193820ca30c1cef7

Code I have:

// Point for filterbound
var roi = ee.Geometry.Point([92.1339935642684,21.19472665398992])

// Center the display to ROI.
 Map.centerObject(roi, 11);
// Function to get NDVI.
var getNDVI = function(image){
  var NDVI = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
  var clouds = image.select('QA60').gte(1024).unmask().not().rename('CloudMask')
  return image.addBands(NDVI).select('NDVI');
  //.addBands(clouds).mask(clouds)
  };

// Gets the current time.
var eeNow = ee.Date(Date.now());

// SET OBSERVATION PERIOD
var startDate = '2018-02-15';
var endDate = eeNow;

// Gets image collection of Observation Period
var series = ee.ImageCollection('COPERNICUS/S2').filterDate(startDate, endDate)
          .filterBounds(roi)
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))

// APPLY NDVI FUNCTION TO OBSERVATIONS
series = series.map(getNDVI)
 Map.addLayer(series.first(), {bands: ['NDVI'], min:0.2, max: 1, palette: ['red', 'orange', 'yellow', 'green']}, 'first')


texts: Getting multiple column values at the same time in QGIS using Python

texts:
I need to filter multiple column values using python in QGIS. I am facing issue while writing expression using QgisExpression(qgis.core).
If one of the column values is NULL then also I need to get the row. Is there a better solution to handle this or any other package to handle the expressions in QGIS?
You can see below an example to illustrate
QgsExpression(f"input='{abc}' OR value='{def}'") # here abc = 1 and def = NULL


texts: I'm attempting to create a simple process with CartoDB, I have a layer with number of polygons, I want to to run a SQL with a multipoint line to report which polygons were intersected.


texts: In a line layer, which is rule-based rendered :  As I can see it, it is only possible to ad simple lines.  I want to ad a marker line.  Is that not possible ?


texts: How to Speed Up ArcGIS 10.1 for Desktop Generating Near Table on Polygon Feature Classes?

texts:
I am trying to create a near table to link two datasets. The first consists of the outlines of about 1.8 million buildings, the second a 3 meter buffer around these buildings, which as been merged together where multiple buffers overlap.

I want to identify the buffer closest to each building. I used the Generate Near Table tool which appears to be working, but is running very slowly. I left it running for 47 hours and it only managed to produce 4700 records. As I have 1.8 million to do I want to know if there is anything that I can do to make this tool run more efficiently.

I've got the 64 bit geo-processing installed and enabled multi-core processing.

I don't think it is the computer which is 6 core @ 3.2 GHz, Windows 7 64 bit, 16 GB of ram, with Solid State Disk.


texts: Density nodes / vertices of polyline by percentage in QGis

texts:
I have a polyline feature class with short and very long features. The goal is to density the vertices of the line features by percentage (of length or original amount of nodes or whatever).
 "curvy" lines should get more nodes by choosing densification by 80% then straight lines with few nodes.
Is there something in qgis?
I just know densification by a adding a specific number of nodes or by certain distance.


texts:
I'm trying to query specific amenities in a specific region on a postGIS database.

I'd like to be able to select multiple amenities in the query, and output their OSM ID, type, lat, lon, & name.

For example:

a nice query on amenity:

SELECT osm_id,ST_X(ST_TRANSFORM(way,4674)) AS X, ST_Y(ST_TRANSFORM(way,4674)) as Y,name
from planet_osm_point where amenity='school' limit 5;

I can't get any results from a bounding box search, and that's part of the problem:

/*with "bbox() ~ way"*/
SELECT osm_id,ST_X(ST_TRANSFORM(way,4674)) AS X, ST_Y(ST_TRANSFORM(way,4674)) as Y,name
from planet_osm_point where ST_MakeEnvelope(-125.0011, 24.9493, -66.9326, 49.5904, 4326) ~ way limit 5;

/*with "way @ st_MakeEnvelope()"*/
SELECT osm_id,ST_X(ST_TRANSFORM(way,4674)) AS X, ST_Y(ST_TRANSFORM(way,4674)) as Y,name
from planet_osm_point where way @ ST_MakeEnvelope(-125.0011, 24.9493, -66.9326, 49.5904, 4326) limit 5;

Here's a gist with a bunch of example queries I've tried to make work on the DBs.


texts: Clip raster with postgis outputting result in clipping geometry extent

texts:
I want postgis st_clip to output clipped result in extent of initial clipping geometry.
Currently, it's crop option can be set to output either in initial raster extent or in extent of intersection between initial raster and clipping geometry.
Here is a test case:
WITH blade AS (
  SELECT
    st_tileenvelope(11,1238, 638) geom
),

butter AS (
SELECT
  st_addband(
  st_makeemptyraster(
    256,
    256,
    st_xmin(t),
    st_ymax(t),
    (st_xmax(t) - st_xmin(t)) / 256,
    - (st_ymax(t) - st_ymin(t)) / 256,
    0,
    0,
    3857),
    '8BUI'::text,
    1) AS rast
FROM
  blade,
  st_translate(geom, 1000, 1000) t
)

SELECT
  st_width(clip),
  st_height(clip),
  st_envelope(butter.rast) AS raster_envelope,
  blade.geom AS blade,
  st_envelope(clip) AS clipped
FROM
  butter,
  blade,
  st_clip(butter.rast, blade.geom, TRUE) clip;

I'm using a random tile as a clipping geometry, then making a test raster by just shifting it by 1000 meters north and east, then clipping it.
Resulting raster size is 243 x 243 and is an intersection between clipping geometry and initial raster.
Image below shows the operation:
Red square is a clipping geometry, green square is initial raster extent, and hashed area is clipping result extent.

I'd like a result to be 256 x 256 pixels and it's extent to look like this:

I've tried to recompile postgis after changing this code in rtpg_mapalgebra.c file a bit, but with no luck yet.


texts:
I found contradictory sources for what's considered in geoid and can't decide on either side because the sources seem about equally reliable. In particular, the question I'm asking is whether acceleration or acceleration due to earth's gravitational force is the same for every point on the surface geoid defines (centrifugal force being considered in the first case but not in the second one).

Furthermore: Let's consider a position http://quicklatex.com/cache3/b2/ql_4f7c7ba7a94c1766ba4e245f67c10db2_l3.png in spherical coordinates. Is there a functional relation http://quicklatex.com/cache3/df/ql_fc5885a7200cf4c05395096db67d04df_l3.png mapping a position in spherical coordinates without the radius to the radius? There should be for geoid to make sense but given earth not being made of a homogeneous material, I can't really imagine there to be a perfect mapping just considering either acceleration or acceleration due to earth's gravitational force. Are the values of http://quicklatex.com/cache3/5b/ql_dbe323ab303f7a84fac2f523dad32f5b_l3.png required to be continuous? Because then it's obvious it works.


texts: How can I get the coordinates of polygons drawn in Google Maps-My Maps?

texts:
I have a custom Google Maps using My Maps which has some polygons (around 18) and I would like to get the coordinates of each polygon without having to place markers on every corner.
This is an example of the polygons I want to get the coordinates. I only want the coordinates of the vertices or the corners, not what is inside.
Is this possible? How I can achieve this?


texts: Turn by Turn Navigation for Google Map on Laptop

texts: How can I generate voice navigation for turn-by-turn directions for Google Maps on my laptop?


texts: How to find which address are enclosing a block in arcgis

texts:
I have two separates shape files. The first one contains blocks and the second one contains roads with addresses. Which is the suitable method in order to find each block from which address is enclosed?

The issue is that I have to assign to every block all the road addresses that are surround every block, meaning that let's say: Block 1 is surrounding by a street, by b street, by c street, by d street, Block 2 is surrounding by a street, by b street, by c street, by E street etc. So I want to create a matrix that it will have all this information...


texts: QGIS Model Designer - Use name of selected Polygon in Output layer

texts:
I am building a simple model using QGIS Model Designer.
Coming from ESRI, I am having a little trouble.
I am clipping buildings of a city by overlaying and iterating over the districts.
The name of the various output files should be like "buildings_[districtName]".
So I need to get the name of the current overlying district from attribute table.
So far, I wasn't succesful. Do you have any advise, how to do it?


texts: Cannot create route from one point to another ArcGIS Android

texts:
i've been trying very hard to create a route between two points(startPoint, endPoint). But i am getting the following error:

Location "Location 1" in "Stops" is unlocated.  Location "Location 2" in "Stops" is unlocated.  Need at least 2 valid stops.  "Stops" does not contain valid input for any route.

Here is my code:

private final String routeTaskURL = "http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Network/USA/NAServer/Route";

@Override
protected void onCreate(Bundle savedInstanceState) {
    super.onCreate(savedInstanceState);
    setContentView(R.layout.activity_main);
    mMapView = (MapView) findViewById(R.id.map);
    mMapView.enableWrapAround(true);
    new Thread(new Runnable() {
        @Override
        public void run() {
            try {
                UserCredentials userCredentials = new UserCredentials();
                userCredentials.setUserToken(token, clientID);
                RouteTask routeTask = RouteTask.createOnlineRouteTask(routeTaskURL, userCredentials);
                RouteParameters routeParameters = routeTask.retrieveDefaultRouteTaskParameters();
                NAFeaturesAsFeature naFeatures = new NAFeaturesAsFeature();
                SimpleMarkerSymbol sms = new SimpleMarkerSymbol(
                        Color.RED, 5, SimpleMarkerSymbol.STYLE.CIRCLE);

                Point startPoint = new Point(36.793653, -119.866896);
                Point stopPoint = new Point(36.795488, -119.853345);

                Graphic graphic1 = new Graphic(startPoint, sms);
                Graphic graphic2 = new Graphic(stopPoint, sms);

                StopGraphic startPnt = new StopGraphic(graphic1);
                StopGraphic stopPnt = new StopGraphic(graphic2);

                naFeatures.setFeatures(new Graphic[] {startPnt, stopPnt});
                routeParameters.setStops(naFeatures);

                RouteResult mResults = routeTask.solve(routeParameters);
                System.out.println(mResults.getStops());
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
    }).start();
}

I've searched the internet. Many developers were/are facing this problem. I've tried all the solutions but none of them worked. I got routeTaskURL from the ArcGIS Routing Sample app. The link which is given in the documentation of ArcGIS maps gives me the 403 error if i open it in the browser.

Note: "token" and "clientID" is declared in the first step and they both are taken from the ArcGIS developers console where i registered my application.


texts: Is GIS a tool or a science? There's a lengthy and dense paper on the matter here. Could someone offer a more digestible summary?


texts: Python :ImportError: DLL load failed: The specified module could not be found

texts:
I am working on Python plugins for QGIS. I set my variables as

PATH :=C:\OSGeo4W\apps\qgis;%PATH% PYTHONPATH:=C:\OSGeo4W\apps\qgis\python

But when I try to run my .py file in IDLE I get an error saying:

Traceback (most recent call last):

File "C:\rt_sql_layer_working\DlgQueryBuilder.py", line 30, in from qgis.core import * ImportError: DLL load failed: The specified module could not be found.

What should I set for the PATH variable??


texts:
I’m doing phenological study using a threshold method in IDL. I have a twelve-band image in 2005 (format is in .tiff not ENVI format). Twelve-band means that monthly images (Jan to Dec) have been stacked in ARCGIS.

I know how to set threshold value on a single image. But I don’t know how to set a threshold value on multiband image to determine the date when growing season starts and ends. Start of the season is determined as the day when the NDVI value is return to threshold value in a given year.

So what I’d like to ask:

Can the stacked image be in .tiff format to be imported into IDL via the envi_open_data_file procedure. Or is an ENVI-formatted file required as input?
What kind of procedure and function should I use to extract start and end of season using threshold method on multiband image? The start/end date of season varies by vegetation type. Basically I want one image showing different places are different phenological date. Is there anyone who can help me for this?

Thanks a lot


texts: Accessing QGIS map canvas as a variable for expressions

texts:
Is it possible to access the QGIS Map Canvas as a variable or expression value to be used in the field calculator or labelling engine?

(i.e. so that you can create an expression that selects/interrogates features that are visible/present in a particular view of the canvas in a similar way that the print composer 'filters the legend by map content')

(QGIS 2.18.1 [osgeo4w], win7)

EDIT:
Potential (fictional) use case 1 - There are 50 polygon's in a layer, currently there are 20 polygon's visible in the current view of the map canvas, I would like to style the polygon's blue if there are <= 25 polygon's visible or style them green if there are >= 25 polygon's visible.

Potential (fictional) use case 2 - There are 50 polygon's in a layer, currently there are 10 polygon's visible, 5 that have an attribute value of 'House' in the "Building Type" field/column, and 5 that have the value of 'Garage' in the "Building Type" field/column, visible in the current view of the map canvas. Return: the value equalling the number of polygon's that are garages ... answer = 5.


texts: Making acres or sq km field column display permanently in ArcGIS Desktop?

texts:
I have a polygon and I calcuate acres or square kilometers.

I know how to do this but every time I bring the polygon to ArcMap. I had to go to their column and set it up the numbers to display look like this 1594267 rather than 1,594,267 ?

I know you can manually set up the column's properties and set it as thousand separators but Is there a way around to make it permanently ?

I have screen shot that you can check this out here.

There are two files before and after.


texts: How to prevent commit to a vector layer during editin session in PyQGIS

texts:
I am writing a plugin that listens to editing events, such as the QgsVectorLayer.editingStarted, QgsVectorLayer.geometryChanged signals. So far all is fine but there are also some validation rules that should be applied at the end of the editing session and commit should be prevented if they fail. The commit should not go through and the editing session and edit buffer should stay unaffected.
On the C++ side, I can call QgsVectorLayer::setAllowCommit() from within the beforeCommitChanges() slot but this is not available from Python, as the QGIS documentation specifically says. Is there some workaround that can be used from PyQGIS to get the same effect?
I tried raising exception from beforeCommitChanges() but this still proceeds with the commit. Subclassing QgsVectorLayer and overriding commitChanges() does not help either as I suspect the commit is triggered from C++ and it does not go through Python at all?


texts: How can I trace an ANALYSIS_EXPORT error when building QGIS?

texts:
When building qgis1.8 using QT 4.7.4, I encountered the following error. I am unable to trace it:

src/analysis/raster/qgsrastercalculator.h:31:24: error: variable ‘ANALYSIS_EXPORT QgsRasterCalculatorEntry’ has initializer but incomplete type

Earlier there was an error for CORE_EXPORT too, i have solved it using DEFINES += GUI_EXPORT= CORE_EXPORT=

Any help is greatly appreciated.


texts: Reassigning raster values by looking for nearest pixel in another raster using QGIS?

texts:
I have two non-overlapping rasters, A (red hues) and B (blue hues). I want to change the value of pixels in B via a nearest-neighbor search in A; i.e. a pixel in B with a value of 15 looks for the nearest pixel in A, finds a nearby A pixel with a value of 30, then the pixel in B changes from 15 to 30.

I can do this by converting both rasters to points and using distance matrix and joins as in Joining attributes from nearest point in QGIS?, but it's slow and kind of cumbersome; is there an all-raster way to do this in QGIS?


texts: Problem in importing raster data from PostGIS to Geoserver

texts:
I am using Opengeo Suite 3.0.2, PostGIS 2.0.3 and PostgreSQl 9.2.

I am trying to import raster data from PostGIS to Geoserver but facing an error.

I follow this tutorial link to import the raster data from PostGIS.

http://dev.horizon.opengeo.org/opengeo-docs/geoserver/tutorials/imagemosaic-jdbc/imagemosaic-jdbc_tutorial.html

All steps are going good.

I used FWtools 2.2.4 to make tiles.
gdal_retile works fine.

These are my XML configuration files...

rasterdata.pgraster.xml

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ImageMosaicJDBCConfig [ <!ENTITY mapping PUBLIC "mapping" "mapping.pgraster.xml.inc">
<!ENTITY connect PUBLIC "connect" "connect.pgraster.xml.inc"> ]>
<config version="1.0">
<coverageName name="slp"/>
<coordsys name="EPSG:4326"/>
<scaleop interpolation="1"/> &mapping; &connect;
</config>

connect.pgraster.xml.inc

<connect>
<dstype value="DBCP"/>
<username value="postgres" />
<password value="postgre" />
<jdbcUrl value="jdbc:postgresql://localhost:5432/spatial_db" />
<driverClassName value="org.postgresql.Driver"/>
<maxActive value="10"/>
<maxIdle value="0"/>
</connect>

mapping.pgraster.xml.inc

<spatialExtension name="postgis"/>
<mapping>
<masterTable name="rasterdata_master2" >
  <coverageNameAttribute name="slp"/>
  <maxXAttribute name="maxx"/>
  <maxYAttribute name="maxy"/>
  <minXAttribute name="minx"/>
  <minYAttribute name="miny"/>
  <resXAttribute name="resx"/>
  <resYAttribute name="resy"/>
  <tileTableNameAtribute name="TileTable" />
  <spatialTableNameAtribute name="SpatialTable" />
</masterTable>
<tileTable>
  <blobAttributeName name="rast" />
  <keyAttributeName name="rid" />
</tileTable>
<spatialTable>
  <keyAttributeName name="rid" />
  <geomAttributeName name="geom" />
  <tileMaxXAttribute name="maxX"/>
  <tileMaxYAttribute name="maxY"/>
  <tileMinXAttribute name="minX"/>
  <tileMinYAttribute name="minY"/>
</spatialTable>
</mapping>

After that I used the DDL generation utility which creates 4 Sql Scripts.

I execute the required script which creates the metadata tables.

now I used java utility to import the image data. That also works fine.

But when I try to use the image mosaic jdbc plugin to import the raster data through jdbc database, it gives an error, i.e.  "Could not list layers for this store, an error occurred retrieving them: Argument "value" should not be null."

I am not able to insert the screenshot right now, so can't show you the screenshot of error.

Thanks

Akhil


texts: DEM generation from ASTER L1B imagery using FOSS

texts:
Is there any (photogrammetric or other) FOSS package to support the generation of a DEM from ASTER L1B products (i.e. combining the spectral bands 3N and 3B)? Are there any practical examples or/and tutorials? It would certainly be preferable if the tools in question integrate easily in (workflows implemented in) well known GFOSS GIS/RS tools.

Relevant Information

Information about ASTER and its products specifications available at:

ASTER User's Guide - NASA
ASTER Level-1B Registered Radiance at the Sensor - LP DAAC
ALGORITHM THEORETICAL BASIS DOCUMENT FOR ASTER DIGITAL ELEVATION MODELS (STANDARD PRODUCT AST14) -- references the PCI software.

Related ASTER data products

On Demand Digital Elevation Model, AST14DEM
On Demand Registered Radiance at the Sensor - Orthorectified
AST14OTH
On Demand Digital Elevation Model & Registered Radiance at the Sensor - Orthorectified, AST14DMO

DEM-generation related works

DEM and orthoimage generation from ASTER L1B images for remote areas
Comparison of automated digital elevation model extraction results using along-track ASTER and across-track SPOT stereo images
ASTER DEM performance
Digital elevation model (DEM) generation and accuracy assessment from ASTER stereo data

Various

http://www.ipi.uni-hannover.de/fileadmin/institut/pdf/dowman.pdf


texts:
How can I get the value of a point that falls between 2 points of another layer?

I have a Shapefile of points with values representing the distance from an origin (blue in the example bellow) at equal interval along a polyline and a Shapefile of points that I want to get the values (red in the example bellow).

I'm thinking of using the IDW interpolation, but I'm not quite sure...

EDIT 1 : Made a more accurate example

EDIT 2 : Would a TIN (linear interpolation) and then a TIN to Raster give me a better result?


texts: GeoServer make layer style to show grid lines?

texts:
Is possible to create a GeoServer layer style so that gridlines (black lines in the image) of the gridset would be visible?

Update1: trying to achieve following view (red lines represent black gridlines in the first image) using only layer style options.


texts: Elevation shading of linear feature compared to surrounding ground level in QGIS

texts:
I have a digital elevation model (DEM - geoTIFF) in QGIS 3.2. There is a linear feature (earth bund/wall) in the DEM that I would like to shade/colour based on its elevation above the surrounding ground level on which it is built. I have provided an example output below copied from another project, with the red area representing a gap in the earth bund and the green representing the desired 2m above ground level. I'd like to recreate this example using my DEM.

Shading the whole DEM based on elevation is not an option due to the size and range of elevations present, as well as the varying ground level along the earth bund. I suspect the method involves clipping the relevant bund area  from the DEM and some clever work in the Raster Calculator that looks at adjacent values, e.g. Pixel X from within the bund area is compared to Pixel Y located 10m perpendicular to Pixel X, outside of the bund area.


texts: I'm working on a project for a conservation plan. Part of the plan requires that areas added for conservation have less than one structure per five acres (these are pretty large, mostly undeveloped areas). I've done kernel density and and point density analyses before but always with many, many more points. In the past I've had problems if there are not enough points. Should I even be looking at density analysis options? Would Thiessen Polygons or a cost surface be a better option?


texts: How do you connect to a SDE instance on EC2?

texts: We are exploring using Amazon EC2 with ArcGIS Server - and none of us has any experience with it. The question "How can we, from within our intranet, connect to the remote SDE instance in the EC2 cloud?" has come up. Do you just do a direct connect and instead of the server name, use the IP of the EC2 instance? Something like sde:sqlserver:<ip_address> instead of sde:sqlserver:HOU-SQLPRD10? Surely connecting to the remote instance is possible, huh?


texts: Help with Meters at Lat calculation accuracy, two methods

texts:
So I'm coding a program (c++), and thing leads to thing, and I need to find out how many meters to a degree of longitude at a specific latitude to a very high precision. Optimally < 1 meter at all latitudes(The program should be able to use this to plot many points given in lat and lon that are sub meter apart.)  pretty much anywhere where there would normally be people (this means everywhere but within a few degrees of the poles).

When researching I came across some formulas to do this. The one that was most commonly available to me produced distances that were ~50 meters off for every 15 degrees of lat compared to what Wikipedia (and other online calcs) said they should be. The second formula seems to be more accurate, but I have no idea WHY its more accurate then the first one.

code/formulas:

double MetersAtLat(double lat)
{
    long double pi = 3.14159265358979323846;
    long double a =  6378137.0; //equatorial radius, meters exactly
    long double b = 6356752.3142;// polar radius, meters exactly
    long double eccentricitySquared = ((a*a) - (b*b))/(a*a);//0.00669437999014?;
    long double rlat = (lat) * (pi/180);

    //method 1; off by ~ 50m for every 15 degrees away from equator. (for some reason)
    long double answ = (pi * a * cos(rlat))/180 * sqrt (1 - eccentricitySquared * (sin(rlat)*sin(rlat)));

    //method 2; off by ~0.5m? at the equator and gets more accurate the more northern lats given?.
    long double tanPsi =  (b * tan(rlat))/a;
    long double psi = atan(tanPsi);
    long double answ2 = (pi * a * cos(psi))/180;

    return answ2;
}

Is the problem with method 1 that it really is that inaccurate or is it really an issue with my computer/compiler/language? and is the second method as accurate as it looks? (are the numbers it returns good?)

edit: I've been using a table on Wikipedia to verify the numbers. (https://en.wikipedia.org/wiki/Latitude)


texts: I cannot find XY tools plugin in QGIS repository.


texts:
How I can use the leaflet popups events through CartoDB. I have seen that some events are working like 'click', 'movestart',but 'popupclose' and 'popupopen' are not working. My goal is to do some tasks when the user closes the popup.

Can be the problem that I am using SVG in the infowindow?

I am using cartodb.createVis or should I use cartodb.createLayer instead?  am I not noticing something?

cartodb.createVis('map', url)
        .done(function(vis, layers) {
            layers[1].setInteraction(true);
            layers[1].on('featureClick', function(e, latlng, pos, data, layerNumber) {
                cartodb.log.log(e, latlng, pos, data, layerNumber);
            });

            var map = vis.getNativeMap();

            map.on('click', function(e) {
            });
            map.on('popupclose', function(e) {
                var marker = e.popup._source;
            });
            map.on('popupopen', function(e) {
                var marker = e.popup._source;
            });
        });
      }


texts: Invalid geometry constructed by EE Python API

texts:
I have a problem with Python API of Google Earth Engine.

I'm trying to create a vector mask of a small area based on NDVI. I take Sentinel-2 data as two bands and pass a rectangle (defined as series of hand-picked coordinates) to constrict the ROI. When I try to run the attached code I get an error message:

  File "C:/files/test.py", line 24, in ndROI
    ndviVectors = ndviMask.reduceToVectors(ROIgeometry, 10, 'label') # At this point the "EEException" is thrown
...

  File "C:\Users\glowa\Miniconda2\lib\site-packages\ee\geometry.py", line 83, in __init__
    raise ee_exception.EEException('Invalid GeoJSON geometry.')

EEException: Invalid GeoJSON geometry.

Apparently the constructor returns an invalid polygon, although when printed it looks fine. How can I make sure, that geometries are correct? What I tried and had no success with were:

different CRS and different parameter setup
different coordinate declarations (3 levels of brackets and two)
different places
even replacing Polygon function with Rectangle function (and changing coordinate format accordingly)

all ending up with the same exception

I'm using GEE Python package installed in Miniconda on Windows 10. Below the problematic code.

import ee

ee.Initialize()

# Define function, which produces a NDVI-based vector ROI
def ndROI(imgName,rectROI,diffBand1, diffBand2):

    # Production of rectangle limiting ROI - it's a ee.geometry defined from coordinates and staying in the GEE.
    ROIgeometry = ee.Geometry.Polygon(rectROI, False, 0, True)
    print ROIgeometry
    # Land looks similar to ice - do NDVI from Sentinel 2 to mask away. Get S2 from the latest time, compute ndvi, export. Will do once.
    s1forNDVI = ee.Image(imgName) # Choose nicest scene by hand
    ndvi = s1forNDVI.normalizedDifference([diffBand1, diffBand2]) # select proper bands of sentinel
    ndviMask = ndvi.gte(0.2).updateMask(ndvi.lte(0.2)) # Mask greater-than-0.2, that is plants have 1, water and ice have 0.

    # Vectorise NDVI
    ndviVectors = ndviMask.reduceToVectors(ROIgeometry, 10, 'label') # At this point the "EEException" is thrown
    ndviResultGeometry = ee.FeatureCollection(ndviVectors).geometry(); # convert NDVI product to feature collection and take its geometry
    return ndviResultGeometry

# Define inputs
imgName = 'COPERNICUS/S2/20150817T095016_20160923T153347_T34UDE'
rectROI = [[20.1149, 53.8063], [20.1369, 53.8063], [20.1369, 53.7913], [20.1149, 53.7913], [20.1149, 53.8063]]
diffBand1 = 'B8'
diffBand2 = 'B4'

# Run function
eeNDVIpoly = ndROI(imgName,rectROI,diffBand1, diffBand2)


texts: Combining two reclassified raster images in ArcGIS Desktop?

texts:
I have two images one is a subset of the other and have reclassified both images in arcmap.

The subset has been reclassified to assign zero value to the pixels while the main image has been reclassified to various pixel values, however the area of the subset within the main image is wrongly classified.

To address this problem I want to merge the two reclassified images so the subset area is assigned zero pixel value in the final output.

I have tried various tools in arcmap including combine, raster calculator, and merge but both give outputs which excludes areas outside the subset.

How do i combine these two images so the final output contains both the reclassified subset and the main image?


texts:
I have succeeded to show Openlayers map in my application, when i show css popup(not Openlayers Popup) map has to blank as the rest of the application. But Openlayers map is not getting blank.

What could be the problem?

Instead of Openlayers if i use Google Maps its getting blank when i show popup.

I am not using any css related to Openlayers, just i included this http://openlayers.org/api/OpenLayers.js


texts: How to display mouse position as tooltip in OpenLayers-2?

texts:
I want map coordinates as mouse over effect in OpenLayers. I am using the following code. However it is showing pixel coordinates:

    map.events.register("mousemove", map, function(e) {
      var position = e.map.x + e.xy.y;
      OpenLayers.Util.getElement("tooltip").innerHTML = position
    });


texts: Map is not shown in Print Composer or in PDF (QGIS 2.18.3)

texts:
I already saw the other questions concerning the same problem, but after trying everything described in the answers there, I still have the problem that my saved print compositions do not show the maps. Not only do they not show in the print composer but also in the printed PDF version there's only white.

As I work with bigger compositions with many maps in one file, it would be great if there is a way to "repair" these.

It was mentioned in Map does not appear upon opening print composer in QGIS that this could happen, if you delete or modify the map layer -
 because it didn't work properly I had to remove my basemap content .

Could this be the reason for this problem and if so, is there a way to repair it?


texts: Gdal: clipping a raster with another raster

texts:
I'm writing a simple utility to crop batches of multi-band geotiff raster files to the same (smaller) area. Using gdalwarp, I can easily crop a file using a single-polygon clipping shapefile:

gdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif

However, the actual area I want to clip to will always be initially defined by another geotiff raster file, not a shapefile. It would be nice if I could use the extent of that raster as the clipping file, but I'm not sure how to do this. Unsurprisingly, the following doesn't work (it doesn't raise an error, it just doesn't produce anything):

gdalwarp -cutline clipper.tif-crop_to_cutline input.tif output.tif

So, my question is, is there a way to supply a raster to gdalwarp -cutline? Alternately, is there another gdal function that can clip a raster using another raster? If neither of these are possible, is there a very simple way to produce a shapefile with a single polygon defined by the extent of a raster?

This code will be wrapped in a more extensive python script, so I can use either command line gdal utilities or any of the python bindings for gdal.

As a side note, I know that I could easily just make a clipping shapefile that covers the extent of my raster in QGIS. I may wind up doing that if I don't find a straightforward solution, but I will ultimately wind up using this utility on dozens if not hundreds of areas as part of a large automated analysis, so I'd prefer not to have a tedious manual step even if it is very easy.


texts: How are the symbols created for building footprints to have front door (small angle seen in first picture)? I have attached what I have thus far - want to convert the bottom picture to the top for the footprints that are very close together?


texts: performance problem with PyQGIS: manipulating attributes

texts:
I have selected a vector layer with attributes slope and area. I want to combine all the features with slope 0 to 5, 5 to 10, ... , 85 to 90. The original features are to be discarded, and new features for each slope-category are to be created (consisting of combined geometry, and combined area attribute).

I am pretty new to python programming, and while the script I created below seems to do what I expect (?), it has performance issues. It begins reasonable fast, but then slows down to a painfully slow speed.

I wonder what the reasen is for this behaviour, and what I could do better (and why) ?

import math
layer = iface.activeLayer()
for slope_cat_nr in range(1, 19):
    print "Slope Category: ", slope_cat_nr*5
    slope_max = 5*slope_cat_nr
    slope_min = slope_max - 5
    if slope_min == 0:
        slope_min = -0.001
    newGeom = QgsGeometry.fromWkt('GEOMETRYCOLLECTION()')
    anzF = layer.featureCount()
    onePercentAnz = math.ceil(anzF / 100)
    print "Nr. of Features: ", anzF
    progress = 0
    fIds = []
    for f in layer.getFeatures():
        if (progress % onePercentAnz) == 0:
            print (100*progress / anzF),
        progress += 1
        if f['slope'] > slope_min and f['slope'] <= slope_max:
            newGeom = newGeom.combine(f.geometry())
            fIds.append(f.id())
    newFeature = QgsFeature(layer.pendingFields())
    newFeature.setAttribute('slope', slope_max)
    newFeature.setAttribute('area', newGeom.area())
    print "New Area: ", newGeom.area()
    newFeature.setGeometry(newGeom);
    layer.dataProvider().addFeatures([newFeature])
    layer.dataProvider().deleteFeatures(fIds)


texts: How to read multiple GRIB files in R and take mean of those rasters?

texts:
I have downloaded multiple monthly mean GRIB files from https://cds.climate.copernicus.eu/cdsapp#!/dataset/ecv-for-climate-change?tab=form for Surface Air temperature , from 1979 to 1990. 144 GRIB files have been saved in my working directory. I know how to read one GRIB FILE (this is of 1979 JAN):
#setting working directory
setwd("C:/Users/Halima Usman/Documents/ICIMOD")
#reading grib files
    pre1=raster('1month_mean_Global_ea_2t_197901_v02.grib',package="raster")
    test=brick(pre1)
    raster1=stack(test)
    plot(raster1) #january 1979

Now, I would like to know how to read all the 144 files in one go and create a single stack of all those GRIB files and take a mean.


texts: How to get a contourline shape into ArcGIS or QGIS from an AutoCAD file

texts:
My question this time is about a request. I need to create an elevation profile and I know I can do that in QGIS 2.4.0 or ArcGIS 10.2. But first I need elevation data, information that is in AutoCAD extension  (.dwg, I presume) and when I use data-> export data and then save it as .shp I don't get elevation information. So I need to extrack the contourlines only or at least with elevation data per line.

Not sure if it's clear but I'll answer anything related to my question.


texts: JSON to Features geoprocessing tool fails when creating SQL Server feature class

texts:
For some reason the "JSON to Features" geoprocessing tool thinks the schema name is the table name and therefore always fails. Seems like it may be a software bug.  Anyone else have this problem?

ERROR 000206: Cannot create FeatureClass 'Database Connections\AP_dbo@server.sde\AP.DBO.Networks_json'. Table already exists [42S01:[Microsoft][ODBC Driver 11 for SQL Server][SQL Server]There is already an object named 'DBO' in the database.] [AP.AP.DBO]

Where "AP" is the DB name and "DBO" is the schema. It should be naming the feature class: "AP.DBO.Networks_json"

SQL Server 2012
ArcGIS 10.4.1


texts: Getting "Error in par(pin = p)" when using spTransform in R

texts:
I am trying to apply a projection to the following dataset from the USGS: https://www.sciencebase.gov/catalog/item/581d052be4b08da350d524ce. After reading in the dataset, the "proj4string" slot shows the data has these CRS arguments:

  +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0

Here is a simplified piece of code that results in the error in question:

library(rgdal)
library(maps)
Highway_longlat <- readOGR(dsn = "SourceFolder")
Highway_merc <- spTransform(x = Highway_longlat, CRSobj = CRS(projargs = "+proj=merc"))
map(database = Highway_merc)

I get this error:

  Error in par(pin = p) :
    invalid value specified for graphical parameter "pin"

I was not able to find any information on this error that seemed related to my situation.


texts: MapServer: using a CRS for which there is no EPSG code

texts:
I am developing a MapServer map file for a map created with Goode's Homolosine projection. The PROJ string is rather simple:

+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs

The IOGP (former EPSG) never assigned a code to this CRS, meaning that it falls into the special case identified in the MapServer manual:

  However, it is often impossible to find an EPSG code to match the
  projection of your data. In those cases, the “wms_srs” metadata is
  used to list one or more EPSG codes that the data can be served in,
  and the PROJECTION object contains the real PROJ4 definition of the
  data’s projection.

Therefore the map file contains the following segments:

WEB
  METADATA
    "ows_enable_request"    "*"
    "wms_srs"               "EPSG:4269 EPSG:4326"
  END # METADATA
END # WEB

# [...]

LAYER
    NAME 'coutries'
    TYPE RASTER
    STATUS DEFAULT
    DATA Countries_Goode.tif
    PROJECTION
        "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs"
    END
END # LAYER

So far so good, however, the GetCapabilities request only lists the CRSs declared in the METADATA section:

<Layer>
    <Name>sample</Name>
    <Title>sample</Title>
    <Abstract>sample</Abstract>
    <SRS>EPSG:4269</SRS>
    <SRS>EPSG:4326</SRS>
    <LatLonBoundingBox minx="-1" miny="-1" maxx="-1" maxy="-1" />

And from this I am not able to issue a GetMap request with the original CRS. What is needed to make the original CRS available to the GetMap request (and listed in the GetCapabilities document)?


texts: I wanted to interpolate the trend values (trend of sunshine hour calculated over 20 years) throughout the study area from the available station points. My confusion is in the selection of interpolation tools I shall use in ArcGIS to do that (IDM, Kriging, Natural Neighbor, spline and all that).


texts: Developing QGIS Action which loads a raster

texts:
I'm using the action code of this question: "How to create a QGIS Action which loads a raster?"
qgis.utils.iface.addRasterLayer('E:\\Plot Sheet Devt\\1974-1984\\1250k\\[% Grid_Ref %].TIF', '[% Grid_Ref %]')

It is very useful for my work and I would need some more functions.
The first one is that the action layer remains selected after the action.
Now after loading raster It doesn't remain selected.
I tried to add this code:
vl = QgsMapLayerRegistry.instance().mapLayersByName('catalogo_CTR')
iface.setActiveLayer(vl)

but It doesn't work.
Other functions could be:

to automatically assign SRC to loaded rasters
to create a group of layer where to load the raster.

UPDATE
Now, way is a linestring.
I tried your code and at the beginning it gave a few hours on the name of the columns that I changed. Then I inserted the centroid of station in the variables of the ST_LineLocatePoint, because it asked me for a point as a second argument. I entered the new code in the application update. everything works but the order is null.
WITH

dumped_ways AS (
SELECT name,
dmp.path1 AS way_part,
dmp.geom
FROM   s_602_ptrc.way,
LATERAL ST_Dump(geom) AS dmp
)
SELECT ROW_NUMBER() OVER(PARTITION BY nway.name, nway.way_part ORDER BY nway._frac) AS id,
station.nome AS station,
nway.name,
nway.way_part,
station.geom
FROM   s_602_ptrc.station
CROSS JOIN LATERAL (
SELECT dumped_ways.name,
dumped_ways.way_part,
ST_LineLocatePoint(dumped_ways.geom, ST_Centroid(station.geom)) AS _frac
FROM   dumped_ways
ORDER BY
dumped_ways.geom <-> ST_Centroid(station.geom)
LIMIT  1
) AS nway
ORDER BY
nway.name,
nway.way_part,
nway._frac
;

I think the problem is in Dumped_ways query. This is the result


texts:
I have some point data which represent daily lat-lon locations of an animal, with an associated timestamp.

I would like to identify all points where STATIONARY = TRUE. A point qualifies as being stationary if a 100km buffer around it overlaps an additional (say) 5 temporally adjacent  points. So if day 10 is my point of interest, I want to ask whether 5 temporally adjacent days are within a 100km buffer of this point. If days 5,6,7,8 & 9; OR days 11,12,13,14 & 15; OR days 8,9,11,12,13 (etc.) are within the buffer, then STATIONARY = TRUE. If, however, days 5,7,9,11 & 13 are within the buffer, but not the alternate (even) days in between, then STATIONARY = FALSE

I think some kind of moving window buffer will provide the solution, but I don't know how to implement this.

I have been trying to get my head around this problem in both ArcGIS and R, but have had no brain waves so far. This is the closest I've got to a solution, but it doesn't quite fit, I don't think: Identification of consecutive points within a specified buffer

Here is some dummy data, which approximates my data structure (though in reality I have twice daily locations (midday and midnight) with some locations missing - but I'll worry about that later)

x<-seq(0,15,length.out=20)
y<-seq(10,-10,length.out=20)
t<-seq(as.POSIXct('2013-07-01'), length.out = 20, by = "days")
data<-data.frame(cbind(x,y,t=as.data.frame.POSIXct(t)))

            x           y          t
1   0.0000000  10.0000000 2013-07-01
2   0.7894737   8.9473684 2013-07-02
3   1.5789474   7.8947368 2013-07-03
4   2.3684211   6.8421053 2013-07-04
5   3.1578947   5.7894737 2013-07-05
6   3.9473684   4.7368421 2013-07-06
7   4.7368421   3.6842105 2013-07-07
... ...         ...       ...


texts:
How can I apply Multiline label option in the label property?

Name field contains a lengthy name. I want to wrap it.

My QGIS version is 1.8.


texts: Creating transect with square polygons using QGIS?

texts:
I would like to create with QGis 2.2 (german version) several square polygons (1km ²) that seamlessly connect to each other, so that a strip ( transect ) of 9km times 1km length results . ( I once attached a picture of how it should look ) - 
http://www.myimg.de/?img=TransektArcMap75b51.jpg

Although I have installed the plugin " Rectangles ovals digitizing ", but can not choose here the edge length . In addition, the following square is to connect seamlessly . In ArcMap the cursor shows me when digitizing contacts with neighboring polygons or their vertexes or lines on , so you could work " seamlessly " .
Furthermore, each of the squares , partially each quadrant , then an attribute (here, an air quality index) to be assigned , in accordance with which the square to be dyed.

Are different colors within a shape based on various values (how it works in ArcMap ) possible?
For information on tools , etc. ..... I use the German user interface.


texts:
I'm trying to develop an offline demonstrator showing hi-resolution bathymetry data using Cesium 1.11 and Geoserver.

I've merged a number of 1/3 arc second areas around Puerto-Rico from NOAA bathymetry viewer using gdal_merge.py, for example San Juan, PR 1/3 arc-second. This appears to setup correctly in GeoServer using layer preview.

I'm loading the data via WCSTerrainProvider, however only half the globe is shown.

var wcsTerrain = new Cesium.WCSTerrainProvider({
    service : "WCS",
    url : "/geoserver/wcs",
    layerName : "puerto-rico",
    waterMask : false
});

I've narrowed the problem down to the WCS "Describe Coverage" response from geoserver.

GEBCO 30-arc second data with global coverage behaves correctly i.e. -180 to +180, -90 to +90
If I force the coverage of the 1/3 second data to be -180 to +180 using gdal_translate then full globe is seen.

My question

Does anyone have any examples of how to use Cesium with a high-resolution subset of terrain data.
Is it possible to use GDAL tooling to "render" the 1/3 data into a GeoTIFF that covers the entire globe without using excessive amounts of disk space?


texts: QGIS: Return if a value/string if a Geometrie had to be simplified by the Simplify Geometry process

texts:
To kind of automate some length and value calculations I have created a  model to just let a Line Layer from an DXF go through it (the layers get selected out).
I had to introduce an extra step to "Simplify" the geometry because some geometries could not be calculated returning a nan value. I have compared the data I get back in QGIS for those geometries compared to what AutoCAD measures, it is very low in difference. With the lowest setting in tolerance 0,000001, it only gives out 0.001 difference in meters in length.
Now considering that it is roughly only one object most of the time, it does not matter much. We always round up, so no problem normally.
But let's assume that it is more than one object, lets say over 10. The difference does stack up.
What can we do to make a value/string be returned in a separate column?
Now I am not really asking to overcompensate the value calculation. If we can get a return of how many geometries that went simplified, that is a big enough step. Then we know how big of the calculation is.
Yes this will be a very rare situation, it will most likely not happen but better to overcompensate than something going wrong.


texts: Generating pivot table in PostgreSQL that ties cities to counties using QGIS (or PostGIS)?

texts:
I have a cities shapefile and a counties shapefile. I have the counties shapefile imported to a PostgreSQL table. I haven't transferred my cities to PostgreSQL yet, but the issue I really want to figure out is how do I create a third table.

Since cities and townships often overlap several counties, I ultimately want a cities table, a counties table, and a city_county pivot table (ID, CITY_ID, COUNTY_ID) so it will play fair with my web application (Laravel Eloquent ORM).

I've been looking through QGIS and can't seem to figure it out, but even if I can accomplish this with SQL in PostgreSQL/PostGIS, that'd be great.


texts: How to logout when there is no cookie found?

texts:
I am using Opengeo Suite and through client SDK I created an application.

When I run the application, browser automatically shows the login panel. After entering the username and password, it takes me to my application. But there is no way to logout and it don't even get a cookie called JSESSIONID. So I didn't get how to make a logout function for this. It logout or what I don't know, but when I close the broswer and enter the url again in new browser, it ask the password again. If I don't close the browser and just close the tab and try to open the application again, it didn't ask the password.

So should I have to ask to close the browser again again to login the application with another user.

Your help will be appreciable.

Thanks

Akhil


texts: Using GDAL to crop an SRTM DEM file from GeoJSON and translating it to bmp

texts:
I have a simple polygon shape in a GeoJSON file.

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              11.67022705078125,
              46.81885778879603
            ],
            [
              12.402191162109375,
              46.81885778879603
            ],
            [
              12.402191162109375,
              47.1514353582905
            ],
            [
              11.67022705078125,
              47.1514353582905
            ],
            [
              11.67022705078125,
              46.81885778879603
            ]
          ]
        ]
      }
    }
  ]
}

I want to use this to crop a DEM file from the SRTM data and export the cropped area into a bmp file that can be later sent to a CNC machine. This particular polygon is in the W020N90 area, so I'm using that data as source.

First, I translate the .geojson to a .shp file:

ogr2ogr -f 'ESRI Shapefile' area.shp area.geojson

And then attempt to use it to crop the DEM file and create a bmp output file:

$ gdalwarp -cutline area.shp -crop_to_cutline -ot Byte -of bmp -overwrite W020N90.DEM out.bmp
Creating output file that is 87P x 39L.
Processing W020N90.DEM [1/1] : 0Using internal nodata values (e.g. -9999) for image W020N90.DEM.
Copying nodata values from source W020N90.DEM to destination out.bmp.
...10...20...30...40...50...60...70...80...90...100 - done.

This results in a "small white rectangle" of a bmp which is not what I want. I would expect a larger bmp that has some visible features, as there are valleys and mountain ranges in the selected polygon.


texts: Adding Shapefile with Symbology from Layer File to ArcMap MXD using ArcPy?

texts:
I have a Shapefile:

C:\DATA\points.shp

A layer file:

C:\SYMBOLOGY\symbology.lyr

And a map document:

C:\MAPS\map.mxd

How can I add points.shp to map.mxd with the properties of symbology.lyr all using python?


texts: Generating a GeoJSON MultiLineString Feature with multiple lines from a GROUP BY clause

texts:
I'm trying to create some GeoJSON directly from PostGIS (v3.1; Postgres 13.2) and have gotten surprisingly far, considering my aversion to raw SQL queries.
Specifically, I want to create a Feature with a MultiLineString. Here's some code:
SELECT jsonb_build_object(
    'type',       'Feature',
    'bbox',       ST_AsGeoJSON(ST_MakeLine(lines.line), 6, 1)::jsonb-> 'bbox',
    'geometry',   ST_AsGeoJSON(ST_Multi(ST_MakeLine(lines.line)), 6)::jsonb,
    'properties', jsonb_build_object(
        'lineProperties',        jsonb_build_object(
            'ids',   array_agg(lines.track_id)
        ),
        'coordinateProperties',  jsonb_build_object(
            'times',   array_agg(lines.time)
        )
    )
)
FROM (
    SELECT track_id, array_agg(extract(epoch from logged_at)) as time, ST_Collect(point::geometry ORDER BY created_at) as line
    FROM coordinates
    WHERE locatable_id = :id
    and user_id = :user_id
    GROUP BY track_id
) lines;

This gives me a result like this:
{
    "bbox": [-165.771653, -82.394628, 150.356346, 85.163403],
    "type": "Feature",
    "geometry": {
        "type": "MultiLineString",
        "coordinates": [
            [
                [-93.644346, 32.831718],
                [-132.701578, -79.916542],
                [74.115335, -20.866482],
                [-79.272251, -27.737648],
                [73.539323, 12.806696],
                [54.706111, -51.355516],
                [-165.771653, -74.537204],
                [150.356346, 85.163403],
                [81.45228, -82.394628],
                [-151.595094, 62.69771]
            ]
        ]
    },
    "properties": {
        "lineProperties": {
            "ids": [
                "93238588-8beb-4174-ac1c-6cd2ac179bdc",
                "93238588-8c67-4dc2-87af-991d5263c02e"
            ]
        },
        "coordinateProperties": {
            "times": [
                [null, null, null, null, null],
                [null, null, null, null, null]
            ]
        }
    }
}

What I am really after, tho is the following:
{
    "bbox": [-165.771653, -82.394628, 150.356346, 85.163403],
    "type": "Feature",
    "geometry": {
        "type": "MultiLineString",
        "coordinates": [
            [
                [-93.644346, 32.831718],
                [-132.701578, -79.916542],
                [74.115335, -20.866482],
                [-79.272251, -27.737648],
                [73.539323, 12.806696]
            ],
            [
                [54.706111, -51.355516],
                [-165.771653, -74.537204],
                [150.356346, 85.163403],
                [81.45228, -82.394628],
                [-151.595094, 62.69771]
            ]
        ]
    },
    "properties": {
        "lineProperties": {
            "ids": [
                "93238588-8beb-4174-ac1c-6cd2ac179bdc",
                "93238588-8c67-4dc2-87af-991d5263c02e"
            ]
        },
        "coordinateProperties": {
            "times": [
                [null, null, null, null, null],
                [null, null, null, null, null]
            ]
        }
    }
}

Basically, I want a MultiLineString with multiple line strings based on their track_id, much like coordinateProperties.


texts: how do I create simple yes/no dialog in qgis python plugin?

texts:
I'm making a QGIS plugin that processes a dataset. I started with pluginbuilder and have a GUI that allows me to set what to do.

As my plugin is making changes, when it detects a certain threshold of changes I want it to pop up a dialog with the old and new values and ask whether to apply changes or skip. How do I do that?

The basic structure so far:

if abs(oldCount - newCount) > 10:
    ????


texts: Saving raster layer to disk using Python: TIFF raster layer loses data type when saved with QgsRasterFileWriter

texts:
I'm trying to save a TIFF image stored as a raster layer. The image is fetched from a server so that its source (the url passed when creating the layer) is a HTTP url: the image should be saved to disk. All of this is done within a plugin in Python. QgsRasterFileWriter does allow me to the save the images, which I've done using code like this:
rlayer = QgsRasterLayer("http://some/path/testimg.tif", "testimg")

pipe = QgsRasterPipe()
pipe.set(rlayer.renderer().clone())
pipe.set(rlayer.dataProvider().clone())

file_writer = QgsRasterFileWriter('C:\\some\\path\\'+'testimg_written.tif')
file_writer.writeRaster(pipe, rlayer.width(), rlayer.height(), rlayer.extent(), rlayer.crs())

This creates a new file on the disk. However, when I open said file in QGIS, they seem visually identical, but core attributes have been modified. The new image has lost the original data type or bit depth (e.g. Float32) is instead saved as 'Byte - Eight bit unsigned integer'. Weirder still, Windows properties still recognizes the image as 32 bit. Original image on the left and processed on the right.

The original images are also single band, but the process seems to force them to be multiband.
My question is this: How do I save a raster layer to disk in Python so that the original attributes stay intact? My goal is simply to save a version that's identical to the one stored as a layer.


texts: Which is better for storage? Geometry column or varchar?

texts:
Objective:

Need to store and retrieve a list of geometries from the database. Basically these are drawn by the user using OpenLayers which I would like to save and present again as a list (kind of like a recent history). When user clicks on a certain item on the list then that feature is drawn on the layer.

My way of going about it:

With my existing code, I would like to construct the feature from the geometry using OpenLayers ie., each of the geometries are sent to the client(as coordinates sets understood by OpenLayers) via AJAX when the client page loads.
I don't want to have a layer and display these via GeoServer.

Question:

With the above in mind, does it matter how I store the geometries in the backend (Postgres with PostGIS)? If yes, which is better, storing the geometries as GEOMETRY column or just as normal strings(ie., varchar)?


texts: How to combine "positive" with "negative" polygons

texts:
I have 2 polygone classes (1 and 2). Inside the polygons are holes (white). Now  I want to combine them with the polygon in witch the holes are. Understand? Sorry for that bad english...

The white holes have to be green or blue, like the sourrounding polygon. Any idea?
Thanks!


texts:
I need a GPS alarm system for my boat, in case of theft. Maybe notifying me any anomaly on iPhone. How are they called "GPS Tracking" ?

I actually don't need to track my boat while sailing, but just knowing if it is moved away from its current mooring spot.

Do I usually need to pay a monthly fee for these services?
How long do the batteries last on such device?

thanks


texts:
Recently, I installed QGIS 3.10.14. I tried to start the application but QGIS would not start after the "setting up GUI" and "Starting Python". I also tried to uninstall QGIS and python from my system and Reinstall them again but still, the problem is not solved.
How can I fix this?


texts: When converting a layer to a KML file, the Name field is not being used for the pop-up title in Google Earth

texts:
I'm using Arc Map 10.7 to create a KML that shows point data representing milepoints along a bunch of roadways. I want the pop-up when you click one of those points to display the milepoint number as the title of the pop-up.
According to ESRI documentation, the title should default to any field named "Name" as long as there are no labels active on the layer in Arc Map. I added a new field to the layer, named it "Name", made it a text field, and used field calculator to populate it with "Mileage " + [milepoint field]. Then I converted the layer to KML using the layer-to-kml tool. I tested it on a shapefile a while back and it worked fine, but for this file the title is still the [Route_ID] field. There are no extra characters and I even moved it all the way to the top of the fields list in properties.
Does anyone know of another reason the KML will not pick up my "name" field as the pop-up title?


texts: Using ArcPy how do you get out a specific category from a shapefile

texts:
I am using ArcPy and wanting to know is there a way to extract a specific category from a shapefile. For example I have a shapefile that I was given that is labled Kansas_pointOfInterest.shp, and I suppose to extract from that only churches (pts I suppose) from that shapefile.

How is this done through ArcPy?


texts: OpenLayers zoom on a collection of points

texts:
I'm new to OpenLayers, but I have a working map with 5 points based on GPS coords working.  however, it's not zooming to the extents of the points.  Wondering what I've done wrong (could very well be a typo, JavaScript is not my forte...).  the "adjust_map_zoom" function is where my problem lies, i believe.
The map is generated, the points are all there and the map is initially set to the mapDefaultZoom setting.  If I change that initialize setting, the zoom changes.  however, as the points spread apart and then come back together, the zoom doesn't change to fit their new boundaries.
var map;
var mapLat = <?php echo $xcenter; ?>;
var mapLng = <?php echo $ycenter; ?>;
var mapDefaultZoom = 12;

function initialize_map() {
  map = new ol.Map({
    target: "map",
    layers: [
        new ol.layer.Tile({
            target: "tile",
            source: new ol.source.OSM({
                  url: "https://a.tile.openstreetmap.org/{z}/{x}/{y}.png"
            })
        })
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat([mapLng, mapLat]),
        zoom: mapDefaultZoom
    })
  });
}

function add_map_point(lat, lng) {
  var vectorLayer = new ol.layer.Vector({
    target: "points",
    source:new ol.source.Vector({
      features: [new ol.Feature({
            geometry: new ol.geom.Point(ol.proj.transform([parseFloat(lng), parseFloat(lat)], 'EPSG:4326', 'EPSG:3857')),
        })]
    }),
    style: new ol.style.Style({
      image: new ol.style.Icon({
        anchor: [0.5, 0.5],
        anchorXUnits: "fraction",
        anchorYUnits: "fraction",
        src: "https://upload.wikimedia.org/wikipedia/commons/e/ec/RedDot.svg"
      })
    })
  });

  map.addLayer(vectorLayer);
}

function adjust_map_zoom(x1, y1, x2, y2) {
  map.getView().fit(ol.proj.transformExtent([ x1, y1, x2, y2 ], 'EPSG:COORDS', 'EPSG:MAP'), {size:map.getSize(), maxZoom:18});
}

if i change the adjust_map_zoom code to
map.getView().fit(source.getExtent(), {size:map.getSize(), maxZoom:16});

or even
map.getView().fit(source.getExtent(), map.getSize());

the zoom does not change at all.
Could someone point out my mistake please?  the bulk of the code was copied from online examples.


texts: Display Multiple Symbols for a single point feature by field attribute vale

texts:
I want to create a point feature layer and add points.
Around the point I'd like to add icons/symbols and labels based on data in fields of an attribute table.


texts: Trying to delete specific features within feature class

texts:
I am working with arcmap 10.2.

I have a database that has a few feature classes in it. Within these feature classes, they have a field labeled as Height. Height is a string field that houses numbers (i.e 54, 54 m, 54 meters).

env.workspace = folder + "\\Base_Intermediate_Data.gdb"
temp_fc = arcpy.ListFeatureClasses()

for temp in temp_fc:

    arcpy.MakeFeatureLayer_management(temp, temp + "_rural")

    if temp.endswith("Polygons_Rural"):

        heights_poly = arcpy.da.SearchCursor(temp, ["height", "OBJECTID"])

        for height in heights_poly:
            if height[0] != None:
                s = re.search(" ", height[0])
                if s:
                    f = height[0].find(" ")
                    if float(height[0][:f]) < 154:
                        arcpy.SelectLayerByAttribute_management(temp + "_rural", "", """ OBJECTID = {0} """.format(height[1]))
                        arcpy.Delete_management(temp + "_rural")

                else:
                    if float(height[0]) < 154:
                        print "OBJ: ", height[1]
                        arcpy.SelectLayerByAttribute_management(temp + "_rural", "", """ OBJECTID = {0} """.format(height[1]))
                        arcpy.Delete_management(temp + "_rural")

I am encountering this error:

ERROR 000732: Layer Name or Table View: Dataset OSM_Polygons_Rural_rural does not exist or is not supported
Failed to execute (SelectLayerByAttribute).

But, I am not sure why. And because of this error, I can't tell if I am correctly deleting the features that meet my criteria.


texts: How can I merge a DEM mosaic with an overlapping bathymetry raster (while keeping only the portions of each that I want)?

texts:
I am tasked with merging a bathymetry raster with a mosaic raster I created from ArcGrid DEM Tiles downloaded from the National Map, for an areal extent spanning the Coastal Lowlands Aquifer System.

The bathymetry raster is needed because the the Mapper only has information for so far off the coast. In order to fill in the Gulf portion of the grid area, the DEM will have to be merged with bathymetry (the raster was created from the GCOOS-STRM30Plus v6 bathymetry dataset).

The bathymetry ALMOST completely aligns with the rectangular extent of the DEM mosaic but not entirely.

Right now the DEM is a rectangle that has values on land, and is completely black as it extends to the ocean.
In contrast the bathymetry raster has values in the ocean (obviously) and is completely white where the land is.

The output I want to achieve is a merged bathymetry/DEM raster, where I keep the DEM values on land and the Bathymetry in the ocean - and somehow blend the contact between them.

I've tried Mosaic to New Raster and so far, no outputs have turned out correctly.
I'm considering using Raster Calculator but I'm not entirely sure what to input to get the results I want. I only need the portion of the bathymetry raster that  portrays the ocean (the actual bathymetry data) and the portion of the DEM mosaic that displays the coastal lowlands.

Is there a way to clip the rest out using the raster calculator and then merge?


texts: QGIS attribute editing doesn't allow values greater than or equal to 10

texts:
I'm trying to modify an attribute field in a shapefile, type double/real by hand.

The existing number, for example, is 100.01 or 123.45, etc. QGIS will not let me type in a number with more than one digit before the decimal place.

I can type multiple digits to the right of the decimal place, but only one to the left.

The current workaround is to copy all the values to a new field, edit those, and then copy them back with the field calculator ( QGIS 3.2.1
)

A possible clue is that the existing field that is troublesome to edit is type: double, type name: real.

The new field is type:

  double, type name: double

Length is 12, and precision is 11. If precision is the number of decimal place to the right, this might partially explain the < 10 issue. But, numbers >= 10 can be stored in the field, I just can't type them by hand, which leads me to believe this is a limitation of qgis, and not the shapefile.

Edit: Tested with arcmap, no problem entering numbers >= 10, so this appears to be a qgis thing. Edited numbers save fine and are visible in qgis. Interestingly arcmap thinks the length and precision are 0 & 0.


texts: Offset between UTM and lat/long coordinates?

texts:
I need to convert some UTM coordinates to lat/long coordinates in order to write those coordinates into some image exifs.

I used the QGIS tool "Add geometry attributes" to do so and it works as intended, but as soon as I create a new project and add both layers together, there's quite some offset between both data sets (see screenshot).

Any guess what's the reason for that? My guess is that the precision of the lat/long coordinates isn't sufficient? But the QGIS tool doesn't give me the option to set the precision. Anything else? Any other way I can do this coordinate transformation in QGIS?


texts:
I am using QGIS with a connection to a PostgreSQL DB with PostGIS. On that DB we have some sharded tables with geography columns. The problem is that I can not see those tables in the browser of QGIS. I do see others tables that are not sharded.

Does anyone have an idea of a possible solution?


texts: Run PostGIS command line functions in pgAdmin?

texts: Is it possible to utilize the PostGIS command line functions, e.g. shp2pgsql, from within pgAdmin III?  I am using pgAdmin to connect to a remote server.


texts:
Through the code I need to calculate the number of pixels of a layer composed by three classes (severe, medium and low), but the code just take into account the medium and low classes.

https://code.earthengine.google.com/90af044b069216a2deb4a2d639b2fee7


texts: Is there a MapCache 1.2.1 (currently last) binary for Windows?

texts:
I use MapServer 6.4.1 installed with OSGeo4W.

I cannot find equivalent MapCache (1.2.1) binary for Windows and I do not know how to compile it. It is supposed to be part of the mapserver suite.

Can anyone point me to a distribution?


texts: Automatic cutting of shapefile based on outline of another shapefile using QGIS?

texts:
I am using QGIS 2.18.

Can anyone explain how to cut the area of an overlapping shapefile?

For example I have a field with a track running through it. The field is one shapefile (a polygon) and the track is another shapefile (polygon) overlaying the field.

How can I cut the track from the field to show the background mapping data, (without manually having to cut around the shape)?

I assume there is a processing tool to do this but I can't get a tool to do what I want.

I essentially want an automatic cutting tool that will cut one shapefile based on another shapefile's outline.


texts: Is there a GDAL function similar to OpenCV's WarpPerspective?

texts:
I't trying to mosaic a GeoTiff raster against another. Initially, they were both georeferenced in WGS-84 CRS.
Using OpenCV, I calculated a homography transformation (a 3x3 projection matrix), that transforms an image to fit the other - I did it using rasters converted to JPEG images. Then, using OpenCV function warpPerspective, I applied my homography to the image and obtained a JPEG mosaic.

Now, I want to apply that transformation matrix to the GeoTiff raster and warp it along with it's geospatial information. However, I can't find an GDAL/osgeo function that would do that: gdal.Warp() seems to be doing something else entirely.

So, my question is: is there a function in GDAL that would accept my raster image and my transformation matrix and output the correctly transformed raster with modified geospatial information?

If not, is there a possible workaround to achieve that?


texts: Did open layers replace the splitWith functionality or just remove it in v3?

texts:
I can't find any documentation on what happened to splitWith.  It used to belong to ol.geom but was only implemented for linestrings.

Does anyone know whether it was removed or replaced and why?


texts: Define Overviews Tool Error when Iterating through Mosaic Datasets

texts:
I have multiple mosaic datasets in a file geodatabase. I need to run the Define Overviews tool on all the mosaic datasets. However, when I try to iterate through the mosaic datasets with ModelBuilder I get an error after the first iteration, meaning it can't move on after the first mosaic dataset.

I had thought the problem was the output workspace since mosaic dataset overviews are are generated and stored in a folder next to where the mosaic dataset resides. Thus the output location is not in Testing.gdb but in a Testing.Overviews folder with sub-folders for the overviews of each mosaic dataset. So, I was attempting to solve the output workspace location problem with this model.

I'm not sure how to solve the output workspace location problem because the it should be
C:\RasterTests\Testing.Overviews\%Name%.Overviews

I am also not sure the output workspace location is even the problem.

ERROR MESSAGE:

  Executing (Iterate Datasets): IterateDatasets C:\RasterTests\Testing.gdb Block2* MOSAIC NOT_RECURSIVE
  Start Time: Thu Jun 29 10:12:40 2017
  Succeeded at Thu Jun 29 10:12:40 2017 (Elapsed Time: 0.51 seconds)
  Executing (Define Overviews): DefineOverviews C:\RasterTests\Testing.gdb\Block270 C:\RasterTests\Testing.Overviews\Block270.Overviews # "196299.9918 6325479.8683 315300.8666 6525992.993" # # 5120 5120 3 NO_FORCE_OVERVIEW_TILES BILINEAR JPEG 80
  Start Time: Thu Jun 29 10:12:41 2017
  2017-06-29T10:12:42.370: Defining overviews
  2017-06-29T10:12:42.371: Computing default overview pixel size
  2017-06-29T10:12:42.443: Defining overviews for pixel size: 1.875
  2017-06-29T10:12:43.057: Defining overviews for pixel size: 5.625
  2017-06-29T10:12:43.452: Defining overviews for pixel size: 16.875
  2017-06-29T10:12:43.484: Defining overviews for pixel size: 50.625
  2017-06-29T10:12:43.511: Defining overviews for pixel size: 151.875
  2017-06-29T10:12:43.539: Defining overviews for pixel size: 455.625
  Succeeded at Thu Jun 29 10:12:43 2017 (Elapsed Time: 2.27 seconds)
  Executing (Iterate Datasets): IterateDatasets C:\RasterTests\Testing.gdb Block2* MOSAIC NOT_RECURSIVE
  Start Time: Thu Jun 29 10:12:44 2017
  Succeeded at Thu Jun 29 10:12:44 2017 (Elapsed Time: 0.45 seconds)
  Executing (Define Overviews): DefineOverviews C:\RasterTests\Testing.gdb\Block280 C:\RasterTests\Testing.Overviews\Block280.Overviews # "196299.9918 6325479.8683 315300.8666 6525992.993" # # 5120 5120 3 NO_FORCE_OVERVIEW_TILES BILINEAR JPEG 80
  Start Time: Thu Jun 29 10:12:45 2017
  2017-06-29T10:12:46.308: Defining overviews
  2017-06-29T10:12:46.308: Computing default overview pixel size
  ERROR 999999: Error executing function.
  Failed to execute (Define Overviews).
  Failed at Thu Jun 29 10:12:46 2017 (Elapsed Time: 0.93 seconds)


texts:
I would like to use attribute value from Workbench as an attribute that says Yes or No to the Fanout options in Workbench settings. So apply value from WB to the advanced fanout settings of WB.

Any idea?

Someone, please?


texts: Redefined topo and shapefile still won't overlay

texts:
I am trying to overlay a topographic map that is a raster and a road shapefile. The road shapefile was originally not defined. I defined the projection to NAD_1983_2011_Contiguous_USA_Albers and I did the same with the raster through its properties window.

I then used the project and made the input the road shapefile and the output coordinate system the raster. This created a new projection file. My assignment says all I have to do is upload the redefined files to ArcMap and they should overlay. They still won't. I checked both properties and the coordinate systems are the same. I don't know how to fix this.


texts: Where can I find a dataset containing locations of all USPS post offices?

texts:
I searched the USPS web site, queried GeoData.gov, did a general Google search and turned up nothing. I tried querying GNIS for feature type post office, but the data is listed as "entered 1993" for most entries, and just looking at some in my general vicinity I can find now-closed post-offices not listed as "historical", currently open post offices listed as "historical", and open post offices with locations listed as unknown. All in all, not a very good approach.

While a national listing would be ideal, a source for just New Mexico would suffice (at least for now.)


texts: Raster Calculator results in NoData values, ArcMap

texts:
I would like to compute the expression displayed in the image with the Raster Calculator tool by using the slope (degrees) derived from the SRTM 30m DEM data. However, it results in NoData values for some of the pixels. The problematic part of the formula that is responsible for the NoData values it seems to be the part surrounded by the red box. This is happening to pixels that have negative sinθ values. In particular, when I solely calculate the sin of the slope raster the output is correct and contains both positive and negative values. However,
when I use the power of 0.8 the negatives sinθ values lead to NoData values.
P.S. None of my DEM and thus slope pixels have NoData values.
I also try to convert my slope raster to unsigned bit and then try the calculations again (without using the Float() function), but still the NoData values were there.
The Power() function does not work as well.


texts: Determine what value used in Editor.OnCreateFeature

texts: How can I determine what Value-Layer (not the layer itself, but the different colors for different attributes) was used in ArcGIS CreateFeature Event, when I insert a new feature using the Editor?


texts: Is it possible to attain a classified raster with several classes in just one QGIS raster calculator operation?

texts:
I'm using qgis' raster calculator and wanted to make a new raster (oraster) from the values of an input raster (iraster). something like:

if iraster>2 then oraster=1; if iraster >4 then oraster=2; if iraster>6 then oraster=3

all in one operation.

Anyone can help?


texts: What am I doing wrong with QGIS v.clean?

texts:
I thought I understood how the GRASS v.clean tool worked in QGIS (in 2.6.1 Brighton) but I keep getting unexpected results and errors with a polygon ShapeFile. I've run the QGIS TopologyChecker with the following rules over the file:

layer must not overlap
layer must not have gaps
layer must not have duplicates
layer must not contain invalid geometries

747 errors were found, all of the gap and overlap type. I tried using the v.clean snap tool to eliminate the overlaps but this happened:

Can anyone point me in the right direction and/or point out what I'm doing wrong? The v.clean tool appears to be removing vertices and actually creating errors. All I want is the slivers in between polygons removed and the polygon holes filled in.

In the end I need a topological ShapeFile to generalize. I've found what I need for generalization here: Generalizing polygon file while maintaining topology in QGIS?


texts: Arcpy Error in adding input field to field map

texts:
I'm doing a script in ArcGIS I'm new to this.  I have to append some field from a table to another, I did this script in Python but I have a problem when it want make the Mapping field with a field of type shape
shape I mean = point, line, polygon ..

The error that I get is this one:

FieldMap: Error in adding input field to field map
Failed to execute (script).

My code is:

# import system modules
import arcpy, os
from arcpy import env

# get parameters
tabellaInput = arcpy.GetParameterAsText(0)

tabellaOutput = arcpy.GetParameterAsText(1)
# Set local variables
fieldMappings = arcpy.FieldMappings()
# add the table to field mapping
fieldMappings.addTable(tabellaOutput)

fieldnames = [f.name for f in arcpy.ListFields(tabellaInput)]
fieldnamesOut = [f.name for f in arcpy.ListFields(tabellaOutput)]

for elem in fieldnames:
    arcpy.AddMessage(str(elem));

if str(elem) not in fieldnamesOut:
    continue

#creat a field map for each field.
fldMap = arcpy.FieldMap()
fldMap.addInputField(tabellaInput,elem) # deve aggiugere questo field

# Set name of new output field
namefield = fldMap.outputField
namefield.name = elem
fldMap.outputField = namefield
# Add output field to field mappings object
fieldMappings.addFieldMap(fldMap)

schemaType = "NO_TEST"
subtype = ""

try:

    # Process: Append the feature classes into   feature class
    arcpy.Append_management(tabellaInput, tabellaOutput, schemaType, fieldMappings, subtype)

except:

    # If an error occurred while running a tool print the messages
    print arcpy.GetMessages()

it does stop inside the

for....

when come a field of type shape

in fact, if I add inside the

 for

the following line of code:

 if elem == 'Shape':
    continue

it works fine... this is pretty strange. Another thing is this, I have more fields in the input table then the target table.  So why in the target don't I see all the fields of the input table?

Should I do a merge and not the append?


texts: Generating tree structure for 3D models?

texts:
I want to create a tool for converting different type of 3D models into Esri's i3s form:

  This repository hosts the specification for Scene Layers which are
  containers for arbitrarily large amounts of geographic data. The
  delivery and persistence model for Scene Layers, referred to as
  Indexed 3d Scene Layer (I3S) and Scene Layer Package (SLPK)
  respectively, are specified.

To begin with, I'd like to build LOD, or layers for a single model, so
how do I calculate a different level of bounding volume for each part of a model?


texts: Calculate difference between values in linear order using ArcGIS

texts:
The best way to describe this without knowing Python or VB jargon is to add these screenshots. I would like to emulate this Excel calculation in ArcMap10.1. MEAS is the distance along my route and ET_Z is the elevation at the route vertex. I just need the difference from the next point down the line from the previous one.

I'm not sure the best approach to this, but figured that there may be simple Python solution for this. Thanks in advance!


texts: Returning error text from QgsProcessingFeedback

texts:
I have a snippet that intersects two map layers (blockLayer  and ogLayer) in QGIS; and prints out when a bad geometry prevents the intersection from occurring. It builds a subclass from QgsProcessingFeedback, based on the answer posted here by @ndawson.
Here's the code:
class MyFeedBack(QgsProcessingFeedback): # create a subclass from the QgsProcessingFeedback class
    def reportError(self, error, fatalError=False): # QgsProcessingFeedback.reportError() is a method
        if "invalid geometry" in error:
            featNum = error.split("(")[1].split(")")[0]
            print("{} has bad geometry".format(featNum))

blockLayer = QgsProject.instance().mapLayersByName("harvestApprovals")[0]
ogLayer = QgsProject.instance().mapLayersByName('Boundary old growth .gdb layer')[0]

#feedback = QgsProcessingFeedback.reportError(self, error, fatalError=False)
#feedback = QgsProcessingFeedback.reportError(error, fatalError=False)

feedback = QgsProcessingFeedback() # instantitate the QgsProcessingFeedback class
print(feedback)

params = {'INPUT':blockLayer,'OVERLAY':ogLayer,'OUTPUT':'TEMPORARY_OUTPUT'}
rslt = processing.run("native:intersection", params, feedback=MyFeedBack())

print("\nDone!")

I'd prefer to create a list to add those numbers to, so as to know which geometries need repair.
I'm just not sure how (or if it's possible) to build a list within an instance of a class, and somehow return that list from the class.
Is there some way I could myList.append(featNum) as the errors are encountered, then get myList back after the process has run?


texts:
I want to query the polygon geometry and get features which are intersecting borders and within the polygon , but i am getting the features within BBOX of polygon .

my query is like this.

service:WFS,

request:GetFeature,

version:1.1.0,

typename:mystore:flightRout,

outputFormat:json,

srsname:EPSG:3857,

cql_filter:INTERSECTS(geometry,POLYGON((12.95269334858834 75.65514701171875
,12.95269334858834 76.11107962890625
,13.455397924014989 76.11107962890625
,13.455397924014989 75.65514701171875
,12.95269334858834 75.65514701171875))
)


texts: Given a raster layer and some coordinates, is it possible to somehow get the color value (in HEX or RGB) of that pixel in QGIS?

texts:
I just want to get the color value based on some coordinates.
Is there a plugin that can do the job?


texts: Calculate the length of geometry line object using python in QGIS

texts:
Calculating the length of a geometry line object using the Field Calculator gives you 2 options:

$length where the length calculated by this function respects both the current project's ellipsoid setting and distance unit settings
length that performs planimetric calculations in the Spatial Reference System (SRS) of this geometry, and the units of the returned length will match the units for the SRS.

Both these functions result different numbers.

That I am trying to do is to make a plugin using python that performs these 2 different functions.

With the following code I can access to geometry length of my selected features and thus the $length in my case:

layer = qgis.utils.iface.activeLayer()
features = layer.selectedFeatures()
for f in features:
    geom = f.geometry()
    print ("Length:", geom.length())

How I can access to length as well?


texts: Canceling stopping edit session in BeforeStopEditing event of ArcObjects?

texts:
I have additional data to be stored when an edit session is saved. If the storing of the data fails (it is located on an internet service) the edit session should not be saved, otherwise data would be inconistant.

Is there a possibility to cancel the stopping of the edit session during the IEditEvents2.BeforeStopEditing event.


texts:
I am currently having trouble implementing a script suggested for a question I previously posted Using arcpy.da.SearchCursor to script snapping one set of points to location of another set of points?
The suggested script is:

arcpy.AddGeometryAttributes_management("POLYGONS", "CENTROID")
arcpy.AddGeometryAttributes_management("POINTS", "POINT_X_Y_Z_M")
# TRANSFER CENTROID XY TO POINTS TABLE
arcpy.AddJoin_management("POINTS", "POINT_ID", "POLYGONS", "PGON_ID")
arcpy.CalculateField_management("POINTS", "centres.POINT_X", "[POLYGONS.CENTROID_X]")
arcpy.CalculateField_management("POINTS", "centres.POINT_Y", "[POLYGONS.CENTROID_Y]")
arcpy.RemoveJoin_management("POINTS", "POLYGONS")
# AND MOVE POINTS
arcpy.CalculateField_management("POINTS", "Shape",
                                "arcpy.Point( !POINT_X!, !POINT_Y!)")

I tried updating the "POINTS' and 'POLYGONS" to my point and polygon shape files. Unfortunately I am unsure what I'm doing wrong/missing in the implementation.

arcpy.AddGeometryAttributes_management("P:\PolygonCentroid.gdb\BuildingFootprints_XY", "CENTROID")
arcpy.AddGeometryAttributes_management("P:\PolygonCentroid.gdb\Situs", "POINT_X_Y_Z_M")
# TRANSFER CENTROID XY TO POINTS TABLE
arcpy.AddJoin_management("P:\PolygonCentroid.gdb\SITUS", "PARCEL", "P:\PolygonCentroid.gdb\BuildingFootprints_XY", "PARCEL")
arcpy.CalculateField_management("P:\PolygonCentroid.gdb\Situs", "POINT_X", "[POLYGONS.CENTROID_X]")
arcpy.CalculateField_management("P:\PolygonCentroid.gdb\Situs", "POINT_Y", "[POLYGONS.CENTROID_Y]")
arcpy.RemoveJoin_management("P:\PolygonCentroid.gdb\Situs", "P:\PolygonCentroid.gdb\BuildingFootprints_XY")
# AND MOVE POINTS
arcpy.CalculateField_management("P:\PolygonCentroid.gdb\SITUS_SNAP", "Shape",
                                "arcpy.Point( !POINT_X!, !POINT_Y!)")


texts: How can I increase Polygon area and draw on center of image in python?

texts:
I have draw Polygon area on image but it will draw on below side corner of image and small area. How can I increase Polygon area and draw on center of image in python.

I am using below code.

import shapefile
from PIL import Image,ImageDraw
r = shapefile.Reader("K:\NEW_GIS\NEW Data\Garpani_NLRMP_Data\Shape files\Garpani_corrected")
shapes = r.shapes()
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 800
iheight = 800
xratio = iwidth/xdist
yratio = iheight/ydist
#print xratio
#print yratio
pixels = []
for x,y in r.shapes()[39].points:
    px = int(iwidth - ((r.bbox[2] - x) * xratio))
    py = int((r.bbox[3] - y) * yratio)
    print px
    print py
    pixels.append((px,py))

img = Image.new("RGB", (iwidth, iheight), "white")
draw = ImageDraw.Draw(img)
draw.polygon(pixels,outline="black", fill="white")
img_new_name="C:/example/Garpani_corrected.png"

img.save(img_new_name)


texts: Talking to a friend of mine in Hawaii and I am from New Mexico. My question to you all is to see if there is a land use/land cover dataset for Hawaii?