texts:
I am trying to generate a time series of degree days, where 30°C is my threshold of interest, using the daily era5 data. The way I calculate the degree days is: DDnonNeg (°C) = Daily Average Temperature °C - 30 °C, negative values are replaced with zero.
The steps I want to do for this are:

Translate the temperature from kelvins to degree celsius
Calculate the daily values of degree days in degree celsius
Generate daily district-specific means
Export a CSV table with the district-specific, daily DD30 for given years.

I think my code does not somehow allow me to use the band DDnonNeg. How can I do this?
This is the link to the district file:
https://code.earthengine.google.com/?asset=users/barborasedova/gadm36_BFA_2
This is is my code:
var dataset = ee.ImageCollection('ECMWF/ERA5/DAILY')
                  .filter(ee.Filter.date('2014-01-01', '2015-12-30'));
var landSurfaceTemperature =dataset.select('minimum_2m_air_temperature', 'maximum_2m_air_temperature');

// Create a geometry representing an export region.
var regions = ee.FeatureCollection('users/barborasedova/gadm36_BFA_2')

//Translate into degree celsius

var tempCelsius = landSurfaceTemperature.map(function(image){
  var props = image.toDictionary(image.propertyNames());
  var Celsius = (image.multiply(0.02).subtract(273.15));
  // Mask where one of min or max temperature has no data
  var CelsiusMasks = Celsius.updateMask(Celsius.select('maximum_2m_air_temperature')).updateMask(Celsius.select('minimum_2m_air_temperature'));
  // to get the mean temperature
  var meanCelsius = CelsiusMasks.reduce('mean').rename('meanTemp');
  // Calculate the DD and make all the negative values 0 (see: https://en.wikipedia.org/wiki/Growing_degree-day)
var DD = meanCelsius.subtract(30).rename('DD');
  var DDnonNeg = DD.where(DD.lt(0), 0).rename('DDnonNeg');
  return ee.Image(Celsius.addBands([meanCelsius, DD, DDnonNeg]).setMulti(props));
});

print(tempCelsius.select('maximum_2m_air_temperature').limit(5))

var proj = tempCelsius.first().projection()

var dailyAverage = tempCelsius.map(function(image){
  var time = ee.Date(image.get('system:time_start')).format('YYYY-MM-dd');
  var districtsMeans = image.reduceRegions({
    'reducer': ee.Reducer.mean(),
    'collection': regions,
        'crs': proj
  }).map(function(feat){
    return ee.Feature(feat).set('Date', time);
  });
  return districtsMeans;
}).flatten();

dailyAverage = dailyAverage.filter(ee.Filter.notNull(['DDnonNeg']));

print(dailyAverage.limit(5));

Export.table.toDrive({
  collection: dailyAverage,
    fileFormat: 'CSV'
});


texts: Dynamic labels between points, along lines

texts:
Goal:
I want to label a set of lines, with the name of the points that come before and after the line segment.

Red labels in the example below indicate what I want to end up with:

What I have:

A point dataset containing all of the points. Each point has a letter
code (one or two letters), and has a line name.
A line dataset, made from the point dataset. As Line Field, in the Points to line tool, I used the lane names. In the example, this is "Line 1" and "Line 2".

Further details:
I want to have labels along each line segment, that refers to its start and end point. In the example, the first segment starts in A and ends in B, and is therefore labeled "A-B". Similarly, in line 2, the segment that starts in BY and ends in BZ is labeled "BY-BZ".

Preferably I would like the labels to repeat when zoomed in, so that they always show which line segment is visible.

Is this possible to solve in ArcGIS 10.4? I have all extensions available, if that matters.


texts: Retrieving Map Composer Items using PyQGIS?

texts:
In a QgsComposition i would like to retrieve the list of all Composer Items.

I don't want to specify item type.

The function should have been composerItems() but the help says it isn't available in Python Bindings (I'm using QGIS 2.18).

Is there a way to access that list ?


texts: Count of the different combination of two attributes in QGIS?

texts:
I've been trying in vain to use the field calculator in QGIS 2.8 to count the combinations of two attributes in two separate columns.

First column is called "Alds_klass" and second is called "bygn_klass". I want to count the different combinations of the two columns, not the sum, but the number of each combination.


texts: "Batch process" in qgis or simple program for resave multiple vector files into files with different CRS

texts:
Right now i,m using QGIS 2.10.1 Pisa, and i wonder if in QGIS there is any feature or plugin to resave multiple vector files that using CRS "A" into vector files using CRS "B" in 1 execution step (Batch process). I mean not save it one by one.

so far i find the closest solution (but no CRS changing) is batch processing for selected features from the toolbox - on pic

But this feature does not have option to change the output CRS. It would be a great help if someone know how to do this in QGIS or in another standalone program that can do this kind of batch processing.


texts: SRID Transformation in PostGIS giving inaccurate transformation

texts:
I am trying to make a transformation in PostGIS from ESPG 4326 to ESPG 5247. The problem is that when the data is imported into ArcMap 10.7 it is clearly in the wrong place (by about 2.5km). I have projected the same layer using ArcMap 10.7 and it overlays the original (ESPG 4326) perfectly. I have noticed both the proj (.prj) files are different (see below). If I copy and paste the ArcMap proj file into the PostGIS version of the same layer it works so it must be something to do with the proj file.
See the code below.
SQL Code Used to make the transformation
ALTER TABLE table_name
ALTER COLUMN geom
TYPE Geometry(POINT, 5247)
USING ST_Transform(geom, 5247);

ArcGIS proj file
PROJCS["GDBD2009_GEORSO",GEOGCS["GCS_GDBD2009",DATUM["D_GDBD2009",
SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Rectified_Skew_Orthomorphic_Natural_Origin"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Scale_Factor",0.99984],PARAMETER["Azimuth",53.31580995],PARAMETER["Longitude_Of_Center",115.0],PARAMETER["Latitude_Of_Center",4.0],PARAMETER["XY_Plane_Rotation",53.13010235415598],UNIT["Meter",1.0]]

Below is the PostGIS proj file (this is the same as the srtext column in the spatial_ref_sys table
PROJCS["GDBD2009 / Brunei BRSO",GEOGCS["GDBD2009",DATUM["Geocentric_Datum_Brunei_Darussalam_2009",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","1056"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","5246"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["latitude_of_center",4],PARAMETER["longitude_of_center",115],PARAMETER["azimuth",53.31580995],PARAMETER["rectified_grid_angle",53.13010236111111],PARAMETER["scale_factor",0.99984],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5247"]]```

---

I went through all the parameters in both projections and they seem to be the same. From what I can tell from looking online, the Hotline Oblique Mercator is the same as the Rectified Skew Orthomorphic Natural Origin projection.

I added an x,y column in the attributes tables and even though they they do not overlay on the map they both calculate the same coordinates. Also when I drop both of them into QGIS, they overlay as expected. So I'm guessing it is something to do with the way ArcMap is reading the prj file that was exported from PostgreSQL. While the parameters are all the same, they seem to be bundled together differently.


texts: Line - Line and Line - Polygon intersection in WGS84

texts:
What is the correct way (which formula, what algorithm, etc) to determine the intersection point of two lines specified in Lat/Lon coordinates on a WGS84 model?

Note that the lines might be very long so a Cartesian approximation is not good enough.

Thanks!


texts: Clipping to Shape using ArcGIS for Desktop?

texts:
I work on a regional level, within that region there are many towns.

Data only exists at a region wide level, I am trying to avoid creating town wide datasets to save time.

I am currently making a series of town wide maps and I've been able to clip to the shape of the town boundary as a workaround to exclude data outside of the town I am working in.

There are occasions when I need label outside the extent of the clipped town bondary, but I when I label outside of the clipped area it isn't shown because I am clipping to town boundary

So, my question is, is there a way to clip to a shape but also show layers and/or label/annotate outside of it of the clipped area?


texts:
I have a Sentinel 2A scene (downloaded from Sentinel AWS storage in UTM projection) and coordinates (a polygon of lat/lng) of a specific plot of land. How can I extract a piece corresponding to that field (preferably in EPSG 3857, so that I could overlay it on Google Maps) a piece corresponding to my plot of interest (or to a given Google Maps tile), using either Python or CLI tools?

I have used gdalwarp to transform the scene image to EPSG:3857, but how can I then transform the lat/lng coordinates of my plot's borders (or the given tile's borders) to the coordinates of the pixels I need to crop? Are there any ready-made utilities or libraries for that?


texts: I'm looking for a vector-grid that has grid-cell boundaries for every tile/zoom-level in the standard x/y/z web-map tiling scheme. Googling has done me no good...anyone know of a source?


texts:
Many debutant GIS users are asking me which one is better QGIS or gvSIG. Since I use Qgis and I don't have any experience in gvSIG, I aways say that both are very nice softwares and the chose is a matter of liking vanilla or chocolate.
But I want to have a better answer for that question, so I would like to hear from people with experience on both of them, what drove you to choose one or the other as your favorite opensource desktop GIS.


texts: How to find the id of the last added feature using pyQGIS?

texts:
Is it possible to programmatically (in python) access the attribute info of the last added feature in QGIS 2.0 before edits are saved? I'm failing to find a way to get the ID and subsequently other attributes of my feature of interest (the one most recently added).

idx = layer.fieldNameIndex('myField')
myList = []
lastEdited = ID THAT I'D LIKE TO RETRIEVE
for feature in layer.getFeatures():
    if feature.id() == lastEdited:
        print feature.attributes()[idx]

From this I'd like to get 'myField' of my most recent addition.

UPDATE: Following lead of this thread, I was able to achieve my goals via the Python console, but am failing to get it to work for my plugin.

def updateTableWidget(self, fid):
    idx = layer.fieldNameIndex('myField')
    myList = []
    for feature in layer.getFeatures():
        if feature.id() == fid:
            myList.append(feature.attributes()[idx])

Via the python console, I was able to connect to my layer like so:

layer.featureAdded.connect(updateTableWidget)

After adding a QgsMessageLog to my function, the attributes were correctly identified. I now have a beginner python issue, however. Calling the function in my plugin as shown below gives me TypeError: updateTableWidget() takes exactly 2 arguments (1 given):

self.updateTableWidget()

The function was called from feature form python file.


texts: In ArcGIS can a feature set use multiple coordinate systems?

texts:
I'm new to learning GIS and ArcMap, but have some confusion about coordinate systems.  I thought I understood that there are two main types of coordinate systems, Projected Coordinate Systems which are in feet or meters, or Geographic Coordinate Systems in lat/long units.

Working on a project in  ArcMap I see that in the Layer Properties, it shows the coordinates systems as a projected system, NAD_83 and also a geographic system, GCS North America 1983 (The block below is a cut and paste from "Layer Properties").

I thought that the coordinate systems were supposed to be mutually exclusive?  Does it make sense that a feature set would be in a projected and geographic system at the same time?? Or does it mean that NAD83 based on GCS North American 1983? Would it be possible to use NAD83 and a differentGCS?

My understanding is that NAD_1983_StatePlane_California_V_FIPS_0405_Feet measures everything in Meters from 34'N and 118'W, so how can it work with GCS that references the Greenwich prime meridian in lat/long?

Projected Coordinate System:  NAD_1983_StatePlane_California_V_FIPS_0405_Feet
Projection:   Lambert_Conformal_Conic
False_Easting:    6561666.66666667
False_Northing:   1640416.66666667
Central_Meridian: -118.00000000
Standard_Parallel_1:  34.03333333
Standard_Parallel_2:  35.46666667
Latitude_Of_Origin:   33.50000000
Linear Unit:  Foot_US

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


texts: ArcGIS API for Silverlight - pass geometries to ESRI service

texts:
How would I test this service to find what it wants?  http://sampleserver4.arcgisonline.com/ArcGIS/rest/services/Elevation/ESRI_Elevation_World/MapServer/exts/ElevationsSOE/ElevationLayers/1/GetElevations

It says geometries.  But does it want a polyline or a point.  I tried just putting a lat/long in the dialog box but no go.

What I have is a result from a densify service

void DensifyService_DensifyCompleted(object sender, GraphicsEventArgs e)
    {
        GraphicsLayer graphicsLayerVertices = MyMap.Layers["VerticesGraphicsLayer"] as GraphicsLayer;
        foreach (Graphic graphic in e.Results)
        {
            ESRI.ArcGIS.Client.Geometry.Polyline p = graphic.Geometry as ESRI.ArcGIS.Client.Geometry.Polyline;

            foreach (ESRI.ArcGIS.Client.Geometry.PointCollection pc in (graphic.Geometry as Polyline).Paths)
            {
                foreach (MapPoint point in pc)
                {
                    Graphic vertice = new Graphic()
                    {
                        Symbol = LayoutRoot.Resources["DefaultMarkerSymbol"] as ESRI.ArcGIS.Client.Symbols.Symbol,
                        Geometry = point
                    };
                    graphicsLayerVertices.Graphics.Add(vertice);
                }
            }

        }
        DensifyButton.IsEnabled = true;
    }

So what I am wanting to do is get each MapPoint geometries from the
    foreach (MapPoint point in pc)

And figure out how to get that into a list or whatever format is needed by the elevation service at the top of the question.  Thanks.


texts:
I have downloaded a 7.5 min DEM . Now while I am trying to mosaic them it gives me an error saying "Input and target dataset should have the same number of bands" . How can I figure out the number of bands for the input dataset.

The DEM I have downloaded is  for Montara Mountain from the following site http://bard.wr.usgs.gov/getDEMSMap.html

I am using ArcGIS 10.0 to  generate a elevation data. So basically using Mosaic to new Raster (Input are two DEMs). Then this output is projected and then using the meters conversion I want to get the elevation data. I tried to get the properties in ArcCatlog but cannot view it. And this DEM data cannot be added in ArcGIS .It throws an error invalid raster dataset.Failed to create a raster layer.


texts: Unable to create group layer inside existing group layer in ArcGIS for Desktop?

texts:
We have recently had ArcMap 10.3.1 installed and have encountered a frustrating bug. When I try to create a new group layer under an existing group layer, the option is greyed out.

I can't find anything to see if it's a setting I need to change, whether it's a bug or maybe the GIS setup in my company. Similarly, the option to paste a layer to another group layer is greyed out. The only way to do it is by pasting via the data frame name and moving the layer or holding Ctrl and moving the layer. I'm pretty sure I was able to do both of the two problems in ArcMap 10.2.1.

Does anyone have any ideas?


texts: PostGIS inserting shapefile with shp2pgsql, wrong coordinate system?

texts:
I'm trying to load a shapefile (link)  into my database with shp2pgsql.

shp2pgsql -D -W UTF8 -I -s 4326 scbel01012011_gen13.shp bel | psql -d stat

no problem so far, but when I'm trying to get some informations from the table by using a coordinate I get no results.
I checked the table, all data are loaded.
Here my query:

SELECT * FROM bel WHERE ST_Intersects(geom,ST_SetSRID(st_makepoint(4.380286, 51.112025), 4326));

Any suggestions?


texts: Can ArcMap 10.2.2 Publish To Server 10.3?

texts: I realize you probably can't but I just want to verify.


texts:
In PostGIS, the ST_length() function returns the length of a linestring in some units, but not in meters.

Is there a function like ST_Length() but in meters?


texts: Create Alpha band in GeoTIFF to make a specific color transparent using GDAL

texts:
I have a simple RGB GeoTIFF with a white background and want to add a 4th band, the alpha, so that if R+G+B=255*3, then ALPHA=0 ; i.e. if all three bands have a value of 255, then make the fourth zero (or some other number).

This question has been asked before on stackexchange but all the answers provided are not working for me.

So far I have tried:

gdalbuildvrt -srcnodata "255 255 255" -addalpha outfile.vrt infile.tif
gdal_translate outfile.vrt outfile-alpha.tif

but this seems to be an OR statement, where if one band is 255, then it is nodata for the whole pixel so a lot of other colors are wiped out as well.

I tried:

gdal_calc.py -A in.tif --A_band=1 -B in.tif --B_band=2 -C in.tif --C_band=3 ^
 --outfile=alpha.tif --calc="255*((A+B+C)<(765))" --overwrite

And this resulted in a single-band raster with 255 in every cell.

I tried the same with --calc="A+B+C" just to see and the all the white cells added up to 253, rather than 765. Is it some kind of modulus thing? No idea...

I have managed to add alpha bands to a white background using QGIS' Raster Calculator using ORs: ((Band@1 <255) OR (Band@2 <255) OR (Band@3 <255))*255 , but I want to be able to do it through GDAL, so that I can do batch processing.


texts: Writing Python script tool to enable ModelBuilder to skip over rasters without data?

texts: This modelbuilder works well until the raster to polygon step. At that step, rasters without features cause an error(error 010151). I would like to put in a step before raster to polygon that causes the model to skip over these "feature-less" rasters. I have very little experience with python and rasters and need some help. I tried to make a row count sequence but it didn't work because the rasters did not have a table to check.


texts: How to refer to a variable from the conditional statement in the ExtractByAttributes function in ArcGIS 10?

texts:
Trying to use the ExtractByAttributes for rasters so that I have multiple attribute values that I want to extract.

This works fine:
myraster = arcpy.sa.ExtractByAttributes("source_raster", "VALUE IN (1,8,10)")

but obviously I can't repeat this procedure for hundreds of layers by manually inserting the desired attribute values for each. Anything I've tried for replacing the numbers with a variable has not worked. It doesn't seem to take a list of strings or integers, nor does it help to strip the brackets[] from around the list, or to have the numbers as a single string. I've even tried with a single number instead of many and I just can't get it to work without inserting the actual number in there. What on earth am I missing? Needless to say I'm a beginner with Python.

Many thanks for any ideas


texts: Auto-Update Color-Ramp to Current Display Extent

texts:
When displaying raster files in Qgis with the color-ramp stretched to the min/max values of the current display extent, the min/max values don't update automatically when zooming or panning to a new area.

The only way i know around it is to manually recalculate the values each time, which can be quite a pain.

I know this functionality to be a basic feature in ArcMap, so maybe I'm just missing something in the display options. So if anyone knows how to get the min/max extent to auto-update, please let me know how to do it.
Thanks!

BTW: I'm currently using Qgis 2.4.0, but I encountered the same issue in all earlier versions I worked with.


texts:
I would like add a feature collection to my map, with a color indicating the group it belongs. Here I show you a screenshot of the results in QGIS I would like to arrive to in google earth engine:

For example, I would like to color the features by country. As an example, imagine I start with the 'states' feature collection:
var table = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");

var states = table.select('ADM0_NAME')

Map.addLayer(states)


texts: Avoid to get rounded corner after a polygon buffer (FME)

texts:
I try to create a buffer for my polygons that leave the angles intact (not rounded).

With postgis/postgresql I can do:
SELECT st_buffer(geom,5,'join=mitre') FROM mypolygon;

With ArcGis I can select the option FLAT in the buffer GUI.

With FME I'm unable to get a correct result. Any IDEA ?


texts: What is naming convention of Tiles generated in GeoServer?

texts:
While the requests to the tiles from the client are in the format {x}/{y}/{z}, the tiles are not named in the same fashion in the file system. The folder names also look different.

What is the naming convention followed by GeoServer while creating the tiles?


texts: Trigger spatial join on insert and update

texts:
I've a problem with a trigger in PostgreSQL.
I can do the query but I don't know to develop a trigger.
The example query is:

select schemaA.table1.value
from schemaB.table2
on
st_contains(schemaB.table2.geom,schemaA.table1.geom)

where geomA is a point geometry and geomB is a polygon geometry.

How can I do this?


texts: How to extract values from metadata in sentinel 2 images

texts:
I am trying to composite an image with the mean value of the METADATA layer "CLOUDY_PIXEL_PERCENTAGE" for every pixel.

My goal is to do this for monthly mean values, but I can't find the way to simply extract the values of this metadata expressed by a percentage (which I know is the format of this layer)

Does anybody know how to?


texts: Converting Grayscale GeoTIFF to PNG using QGIS results in all black image

texts:
I have a GeoTIFF in QGIS that I created by merging some images inside QGIS.
I use Raster>Conversion>Translate. And save it as a PNG.
I then open the image in both, Photoshop and Windows Photo Viewer and it looks entirely black.
Anyway to solve this?
I would like to use the image to print/use it in a document.
Edit: Here is the gdal info:
Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
AREA_OR_POINT=Area
Band 1
STATISTICS_MAXIMUM=1079
STATISTICS_MEAN=252.34346333641
STATISTICS_MINIMUM=1
STATISTICS_STDDEV=118.63833194341
Dimensions
X: 5400 Y: 3600 Bands: 1
Origin
-63,-37.6667
Pixel Size
0.000277778,-0.000277778
No Data Value
0
Data Type
UInt16 - Sixteen bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
-63.0000010703059985,-38.6666649293319793 : -61.5000010691059984,-37.6666649285319792
Band
Band 1
Band No
1
No Stats
No stats collected yet


texts:
I am trying to do a very simple interpolation with soil nitrate data. I wanted to do IDW because points close to each other should be most similar. I used the IDW interpolation function in ArcMap (image 1), and got some strange results.

The layer is stretched/elongated on the vertical axis (image 2). In fact, when I try kriging and other methods from different packages, its still elongated. I'm still new to this and thought I understood the input variables.

Is there a variable or parameter that I am missing which over emphasizes the vertical relationship?

Info on data source:
Geographic Coordinate System:   GCS_WGS_1984, Datum:    D_WGS_1984, Prime, Meridian:    Greenwich, Angular Unit:    Degree

Info on tiff layer: Spatial Reference: Polar_Stereographic, Datum: D_WGS_1984

Edit: When the data are projected using polar stereographic, they no longer correlate to the map and do not represent the original sampling scheme (images 3 and 4)

Image 1:

Image 2:

Image 3:

Image 4:


texts:
I have the following version on my machine. But the MapScript is not installing.

Can any one help?


texts: OpenLayers layer text popup override code not work

texts:
I have some problem with the popup code in OpenLayers, using Text layer (yes i know there already GML layer out there with text format, but i need to modify existing code), the problem was popup not removed when i open new popup (by clicking another marker)

OpenLayers.Feature.prototype.popupClass = OpenLayers.Popup.FramedCloud;

        OpenLayers.Layer.Text.prototype.markerClick = function(evt) {
            var sameMarkerClicked = (this == this.layer.selectedFeature);
            if (this.popup) {
                this.popup.toggle();
            } else {
                this.layer.map.addPopup(this.createPopup(true));
            }
            this.layer.selectedFeature = (!sameMarkerClicked) ? this :
null;
            OpenLayers.Event.stop(evt);
        }

code by Christopher Schmidt : http://www.mail-archive.com/users@openlayers.org/msg03869.html


texts: Setting variables to None in GDAL/Python

texts:
Is there a rule I can follow to know when I need to set something to None at the end?
I had a script where I was looping through a shapefile layer. I had only set the variable with the ogr.Open(file) to None at the end.
However, when I ran the same script again, I realized it was retaining that layer information and different times running the script would give me different results until I also set the layer to None at the end.
I'm not quite understanding what times I need to do this.


texts: Creating geopackage from multiple raster files?

texts:
I've got a bunch of geotiffs' that I would like to combine into a single geopackage. Using QGIS, I've only managed to export one layer as a geopackage. I've tried merging the layers and created a virtual layer but then QGIS messes up the colors.
The raster files are also arranged in a Arcgis Mosaic dataset (file geodatabase) if that is of any help.

Currently, I do the following (which probably isnt the fastest or most efficient way):

Generate a .txt-file containing the filenames of all the geotiffs (dir /b /s *.tif > list.txt)
Use gdal_merge to merge all the geotiffs from list.txt into one large geotiff (gdal_merge -pct --optfile list.txt -o new_huge_geotiff.tif)
Use gdal_translate to convert my new, huge, geotiff into a geopackage (gdal_translate -of GPKG new_huge_geotiff.tif new_tiny_geopackage.gpkg)
use gdaladdo to build pyramids in my new geopackage, which I guess also could be done in the previous step when creating the geopackage. I'm just such a newbie. (gdaladdo -r cubic -oo tile_format=JPEG new_tiny_geopackage.gpkg 2 4 6 8 16 32 64)

The process of merging all small geotiffs into one huge takes a lot of time - I wish I could just add them all to a geopackage at once. I've tried using --optfile and the generated txt-list but gdal_translate doesnt like that command.


texts: R Overlay points and polygons with a certain degree of tolerance

texts:
Also asked here
Using R, I would like to overlay some spatial points and polygons in order to assign to the points some attributes of the geographic regions I have taken into consideration.
What I usually do is to use the command over of the sppackage. My problems is that I'm working with a large number of geo-referenced events that happened all over the globe and in some cases (especially in coastal areas), the longitude and latitude combination falls slightly outside the country/region border.
Here a reproducible example based on in this very good question.
## example data
set.seed(1)
library(raster)
library(rgdal)
library(sp)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(0.30*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts1 <- spsample(p2-p, n=3, type="random")
pts2<- spsample(p, n=10, type="random")
pts<-rbind(pts1, pts2)

## Plot to visualize
plot(pts, pch=16, cex=.5,col="red")
plot(p, col=colorRampPalette(blues9)(12), add=TRUE)

# overlay
pts_index<-over(pts, p)

# result
pts_index

   ID_1       NAME_1 ID_2           NAME_2 AREA
1    NA         <NA> <NA>             <NA>   NA
2    NA         <NA> <NA>             <NA>   NA
3    NA         <NA> <NA>             <NA>   NA
4     1     Diekirch    2         Diekirch  218
5     3   Luxembourg    8         Capellen  185
6     3   Luxembourg    9 Esch-sur-Alzette  251
7     3   Luxembourg    8         Capellen  185
8     1     Diekirch    1         Clervaux  312
9     1     Diekirch    2         Diekirch  218
10    3   Luxembourg    8         Capellen  185
11    1     Diekirch    2         Diekirch  218
12    1     Diekirch    3          Redange  259
13    2 Grevenmacher    7           Remich  129

Is there a way to give to the over function a sort of tolerance in order to capture also the points that are very close to the border?
NOTE:
Following this I could assign to the missing point the nearest polygon, but this is not exactly what I am after.


texts: QGIS python file error for fetching script info

texts: Currently i am working on automated map generation using Python coding in QGIS. i saved my all code but wile fetching the file into Script runner plug in i am getting error that you moved or deleted files. But i saved those all codes. whether it is issue of plug in or my code? I am just started to learning Python coding in QGIS theretofore anyone help is really valuable


texts: Large GeoJSON file in Leaflet performance issue

texts:
I have a large GeoJSON file (now it is about 22MB with 4K polygons but it is going to contains about 10k polygons)

I load the file in Leaflet using L.GeoJSON.AJAX, on each feature I attach two handlers (click and contexmenu) and the puropose of the map si to select polygon read and edit some properties and save the selected polygon as project in a MySQL database.

When I load the map after loading the base polygons from the large file I load every saved project from MySQL database.

For now all works fine and It's quite fast but the GeoJSON file is going to grow so my question is: there is a better way to load the base polygon?

It will be faster If I store the base GeoJSON in the MySQL database?


texts:
I'm a new user of GRASS 6.4.3 and I wanted to import ESRI ASCII grid.

Raster file was unzipped, and in r.in.arc window I specified name for output, storage type etc.

After it has run, raster didn't appear on Display, and in Command Console there was next message:

Command 'd.rast map=39_03@Slovenija_WGS' failed

Can anyone help me?


texts: Data driven pages: merge multiple pdfs based on string per wildcard

texts:
I am trying to figure out how to automatically merge pdf files based on a stringpart (*wildcard) in the filename, without loosing selectable layers within a pdf.

I created 3 series of layout pages with the data driven pages tool. All series include the same index features but there are different layers/raster datasets in the background (topics are: landuse, climate change, site preservation).
For each topic and index feature I have already exported several multiple pdfs which include the topic and index page name.
Now I would like to merge the pdfs for each feature, like:

landuse_feature1.pdf + climatechange_feature1.pdf + sitepreservation_feature1.pdf = map_feature1.pdf

landuse_feature2.pdf + climatechange_feature2.pdf + sitepreservation_feature2.pdf =
map_feature2.pdf

...

landuse_(*wildcard).pdf + climatechange_(*wildcard).pdf +
sitepreservation_(*wildcard).pdf = map_(*wildcard).pdf

After a lot of research I found a free tool in the internet called "pdftk" which seemed to work just fine, but after merging the pdfs all selectable layers/layer informations were gone.

So, I know that there is a way to append single pdfs to a multipage pdf with arcpy, but since I have more than thousands of pdfs to merge, this is not an option for me.

Does maybe anybody had to struggle with this problem before and knows if and how it is possible to use a wildcard for merging pdfs with arcpy?


texts: Reprojecting xy tiles from mercator to Gauß-Krueger in ubuntu terminal?

texts:
I have exported a map from tilemill as mbtiles and converted them to xyz tiles with help of mbutil. The tiles should be used (for learning purpose) with mod_tile or GeoServer and displayed in OpenLayers.

in OpenLayers already exists a bunch of layers using the gauß-krüger projection.

Do you know a program which runs under ubuntu-terminal to reproject these tiles from Mercator to Gauß-Krüger projection to have only one used projection?

Inspecting a raster tile with gdalinfo gives me following info:

Driver: PNG/Portable Network Graphics
Files: /Users/srv/Desktop/tiles/11/586/784.png
Size is 256, 256
Coordinate System is `'
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  256.0)
Upper Right (  256.0,    0.0)
Lower Right (  256.0,  256.0)
Center      (  128.0,  128.0)
Band 1 Block=256x1 Type=Byte, ColorInterp=Palette
Image Structure Metadata:
NBITS=1
Color Table (RGB with 1 entries)
0: 245,253,240,255


texts: Getting writer in FME to output as appointed filetype while using FilenamePartExtractor?

texts:
My goal is to merge all shapefiles in the same folder and outputting them as a gml with name of the folder.

So it looks like this:

Main Dir

Set1

roads.shp
buildings.shp

Set 2

roads.shp
...

I managed to merge the files together on this way:

Using the Shapefile Reader and setting the input to a folder with subdirectories. Then, using the FilenamePartExtractor transformer I extract the name of the folder (set1/set2/..).

Then I write to GML with the name of the folder and fanning out the GML with this name as input. The output is a file without extension, but with the correct name.

Changing the filetype to .gml confirms they have merged succesfully. But why is it not automatically writing as .gml in FME?

nb. it does output a correct .xsd file.

My navigator in FME:

And the output in my \Testdat_GML_Output\ folder:


texts: Are there any examples of WCS GetCoverage requests done using AJAX?

texts:
I am looking through OpenLayers to try to figure out a consistent format for generating views of coverages that can satisfy any WCS server, but mainly I'd like to get GeoServer coverages.

While I've read that OpenLayers does not support WCS, I see in the API that there is OpenLayers.Format.WCSGetCoverage constructor, but I'm not certain how to go about using it.

I'm not sure what to pass in to the function to get the coverage and ultimately I will parse what I retrieve. I know the call below is incorrect, but I'm using it as a mild example. Ultimately based on the configuration of my service I am going to use a DescribeCoverage type of function to retrieve parameters and then subsequently use the GetCoverage upon success of retrieving metadata. Then I need to find a way to parse the XML/GML result, but from what I can tell OpenLayers may not have any parser for that schema.

 $.fn.WcsGetCoverage = function(baseUrl, coverage, onSuccess, onError){
    if(baseUrl.endsWith("wps?")){
       baseUrl = baseUrl.replace("wps?", "wcs?");
       $.ajax({
       type: "GET",
       url : baseUrl,
       data: {
            service : "wcs",
            version : "2.0.0",
            request : "GetCoverage",
            coverage : coverage,
            crs :, //epsg_code,
            bbox :, //minx, miny, maxx, maxy,
            WIDTH : , //output_width: Width in pixels of map picture. One of WIDTH/HEIGHT or RESX/Y is required.
            HEIGHT : , //output_height: Height in pixels of map picture. One of WIDTH/HEIGHT or RESX/Y is required.
            RESXx: , //When requesting a georectified grid coverage, this requests a subset with a specific spatial resolution. One of WIDTH/HEIGHT or RESX/Y is required.
            RESYy: , //When requesting a georectified grid coverage, this requests a subset with a specific spatial resolution. One of WIDTH/HEIGHT or RESX/Y is required.
            FORMAT: , //output_format: Output format of map, as stated in the DescribeCoverage response.
        },
        success : onSuccess,
        error : onError
    });


texts: How to assess geocode service accuracy / recommendations for geocoding services

texts:
I've looked at some similar questions on here, but we have very specific requirements and do not have the technical expertise to figure out how to map other questions to our situation.

We are willing to pay enterprise prices, but we need very high accuracy rooftop geocodes. We would be submitting individual CASS-validated addresses, as many as 15k an hour, and requiring low-latency results. We require a very high uptime. For TOS, we are not a profit-making entity, but we require the rights to store the results. Preference would be given to services that are better able to keep up with new development. Domestic addresses in the State of Wisconsin are required, other domestic addresses and international addresses would be nice. We have no need for reverse-geocoding. We have a use for batch-geocoding but only if expected accuracy is very high (our individual geocodes are human-verified; batches would not be). Responsive customer support is a big plus. Ideally the cost would be less than $5000/year, but accuracy and availability are the two biggest concerns.

I've looked at several companies that provide these kinds of services, and the problem is that our GIS specialist has left, and we lack the technical expertise to evaluate the accuracy of these geocodes. In the interim we are using Bing, but have received many complaints about the accuracy of its geocodes.


texts: need help for script checking topology using ogr in python

texts:
My main goal is to create a new feature indicating if the line features intersect each other as well as if they are not connected to a feature.

My problem is how to identify the intersecting and not connected points and converting them to another geometry layer.

The only thing I came up with is:

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')

inDatasource = driver.Open('test.shp', 0)
inLayer = inDatasource.GetLayer()

#loop through features in the layer
feature = inLayer.GetNextFeature()
while feature:
    #condition
    feature.Destroy()
    feature = inLayer.GetNextFeature()

#check topology

outData = driver.CreateDataSource('TopologyResult.shp')
outLayer = outData.CreateLayer ('TopologyResult', geom_type = ogr.wkbPoint)

fieldDefn = ogr.FieldDefn ('id', ogr.OFTInteger)
outLayer.CreateField(fieldDefn)

featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(point)


texts: Where can I find shapefile (free/open source) for all countries at Admin 2/3 level

texts:
I need to download shapefiles for all country from which I need to generate TopoJSON for around all countries.

Currently I am using data from naturalearthdata.org which only provides Admin levels 0 and 1 for all countries. I need for Admin levels 2 and 3 which covers towns/cities/districts details.

I have tried diva-gis.org I had to download country by country from drop down which is taking lot of time.

Is there any site I can download shapefile for admin 2,3 for all countries as single shapefile or country specific shapefiles. (free/open source)

I will integrate this map/ TopoJSON inside d3.js /Java /j2ee project.


texts: Create Widgets for specific fields of a feature

texts:
I need to create a widget were some fields of a feature are listed and are editable.

I used QgsAttributeForm for this, but this widget is too big and has to many fields.

I need to provide a QCheckBox as well before every field to mark it as selected or deselected.

I took a look at the QgsFeature class and its fields() method. The QgsField class has a type() Method for determing the type of a field, but unfortunately I have not yet found an easy to create a widget from a field type.

I took a look at the QgsFieldWrapper classes, but they don't seem to work in PyQGIS because they are abstract classes.

Is there an easy way to create a default widget from a feature field, e.g.
widget = FieldWidgetWrapper(field)
that I can use to generate these controls?

I know I can create custom UIs for AttributeForms, but this will not work, because I need to be able to have more than one feature in this widget represented by specific fields.

Currently I display everything in LineEdits, but this is not working, especially for DateTime objects. It would be cool if the widgets can be created automatically from a field like in the QgsAttributeForm class...


texts: Extract information from gdalinfo/ogrinfo

texts:
In Bash scripts, I use gdal(srs)info / ogrinfo to get information from raster / vector files.

To use this information in my scripts, I'm extracting information using grep

... e.g. like this (to receive ulx)

gdalinfo raster.tif | grep "Upper Left" | grep -Po '[-0-9\.]+(?=,)'

... or this (to receive EPSG Code)

gdalsrsinfo raster.tif | grep -Po '(?<=AUTHORITY\["EPSG",")\d+(?="\]\]$)'

That works, but I'm not very happy with that solution because it's very dependent on the way gdalinfo works. If anything changes, my scripts might break.

Does anyone know of any better solution to that ?


texts: Dynamic scale bar, coordinates and grid for map frame in Mapinfo map layout

texts:
I am looking for a dynamic scale bar and coordinates (grid) for map frame in the map layout of Mapinfo.

To make it clear I am talking about separate scale bar in layout which is not part of map frame. Of course with possibility to have it automatically connected with map frame so when the change of map is made it also redraw the scale bar and coordinates around map frame (and/or grid).

But I don't see these options anywhere and after some searching it seems it is possible to do it simply and dynamic only with extension tool like Encom Discover.

Am I missing something?

1. Dynamic Graphic scale bar

In MapInfo Professional version 10.0 you can create a dynamic scale bar as part of your map.

Close, but... I don't want it as part of map window, I want it separate. Only thing which come to my mind is make the map frame bigger containing the place where I want scale bar and mask the rest around with white rectangles.

Any better idea for separate dynamic Graphic scale bar?

2. Coordinates and grid for map frame

- no default solution

- free mapbasic tool - from J.Wilcock called Layout Coordinate Grid - very close - allows to draw grid lines and place coordinates around frame. But it is not dynamic, also settings are basic, not many options, maybe having source code, it would be possible to extend it. At least for few more options and possibility to delete all lines and texts at once (again no group function for objects in layout) for redrawing.

Do you know any better way or free tool for coordinates and grid for map frame?


texts: WFS open layers, doesn't return a field result when called

texts:
I've a problem when I call a WFS service, I want to select a file (its attributes) called "DICOFRE".

The service is OK in the institution that produces it, but in my institution when I call the service I have problems in that part of the code (and the code is same as the other institution code)...

Maybe I should change something in the code? If yes what am I doing wrong... if I'm calling the service to my organization and doesn't return the results of the field...

Code:

var wfsProtocol = new OpenLayers.Protocol.WFS.fromWMSLayer(layer_freg);
            var response = wfsProtocol.read({
                filter: pfilter,
                extractAttributes: true,
                infoFormat: "application/vnd.ogc.gml",
                //resultType: "hits",
                callback: function (resp) {
                    if (resp.features == null ||
                        !resp.features.length <= 0) {
                        alert("Não foi possivel gerar DICOFRE!");
                        //layer_editor.destroyFeatures(event.feature);
                        //return false;
                    }

Can somebody help please?


texts: Why does QGIS crash when merging shapefiles?

texts:
This is a repeatable problem. I have created a shapefile using a csv. I am trying to merge it with an existing shapefile, both of which contain only points. (Using the Vector | Data Management | Merge Shapefiles option.) I have double-checked that the shapefile type is correct, and the coordinate reference system matches. As soon as I click OK to accept the result, the message appears "qgis.exe has stopped working" etc., the usual Windows crash message.
Things I have done to reslove this:
1. I have a friend who has the same software and the same project files (we are collaborating on this project). I sent him my original csv file with the new data in it, and he was able to successfully merge the files at his end. So I know that the csv file is OK.
2. Uninstalled and re-installed the program. Identical result, program crashes every time as soon as I go to accept the merging of the shapefiles.

Using Windows 7.

Does anyone know how to resolve this issue? If not, is there another way to bring new data into the project? I don't want to use the "Add Feature" option as it is too difficult to get the co-ordinates correct with just hovering over the map.


texts:
I got the following GeoDataFrame taken from a CSV file and after some slincing and CRS and geometry asignment

    ctf_nom         geometry                                    id
0   Prunus mahaleb  POINT (429125.795043319 4579664.7564311)    2616
1   Betula pendula  POINT (425079.292045901 4585098.09043407)   940
2   Betula pendula  POINT (425088.115045896 4585093.66943407)   940
3   Abelia triflora POINT (429116.661043325 4579685.93743111)   2002
4   Abies alba      POINT (428219.962044021 4587346.66843531)   797

I've converted the geometry from a str trough:

from shapely import wkt

df['geometry'] = df['geometry'].apply(wkt.loads)
df_geo = gpd.GeoDataFrame(df, geometry = 'geometry')

and asigned a crs by :

df_geo.crs = {'init' :'epsg:25831'}
df_geo.crs

when I'm trying to save again the reduced geodataframe by  gdf.to_file() function it returns the following attribute error:

out = r'my_path\df_geo.geojson'
df_geo.to_file(out, driver = 'GeoJSON')

It's even the same when I try with shp extension...

AttributeError: 'Series' object has no attribute 'has_z'

How can I solve this?


texts: What's the difference between the analytic and visual images available through the Planet API?

texts: When I download satellite images from the Planet Labs Scenes Explorer I have the option of an analytic or visual image. What's the difference?


texts: ArcGIS Engine, How to change IMarkerSymbol size in runtime?

texts:
I'm using ArcGIS Engine 10 with C# Winforms to visualize some elements on a globe control.
All elements are of type IMarkerElement which have symbol of type IMarkerSymbol.
Now I want to add two buttons to increase/decrease element's size.
I tried this two methods with no success

1.

((element as IMarkerElement).Symbol as IMarkerSymbol).Size = newSize;

2.

((element as IMarkerElement).Symbol as IMarker3DPlacement).Size = newSize;

I noticed that I can't change IMarkerSymbol's properties values after creation (like size, offset, depth, ..etc.)?
How to change this values?


texts: Gdal and module import error when creating an .exe application with pyinstaller

texts:
I’ve been struggling since three days in a row to create a windows executable application from my python application. I have managed to get results but I’m not at all satisfied with the workarounds I did to make it work…

I am using qgis library and its Python 2.7 interpreter in which I was able to install the freezing tools. I used two different freezing tools : cx_Freeze and pyinstaller and I got the same errors.

To start building executable application with those tools I used a command prompt that I opened with this a bat file (OSGeo4W.bat) situated in « C:\QGIS ». When I have the build folder created with cx_Freeze (or the dist folder with pyinstaller) I always try to launche the executable application with a commend prompt without the previous environnement variable initialisation because I want a standalone application where the user doesn’t need to install qgis library on his computer.

The first error I always get is :

ImportError : No module named QtNetwork

I do not need this module but for the sake of the application, I modified my python code and added the following line :

from PyQt4 import QtNetwork

I will then create a new .exe application and launch it again and I find other missing modules, at the end I changed my python code to :

from PyQt4 import QtNetwork, QtXml, Qsci

Finally I do not have any further missing module but I’m not sure that I realy solved the issue here… ?
gdal folder:

The second error I always get after "resolving" the import module issue is :

ERROR 4 : Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

This seemed to be the same type of error than the previous one so I tried adding this code in my python code and I copied the gdal directory in my .exe application :

gdal_prefix=str(cur_dir)+"\\gdal"
os.environ["GDAL_DATA"] = gdal_prefix

When I start the .exe application I still get the same error :

ERROR 4 : Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

The only way I managed to find a workaround was to create a « Launch_app.bat » file with the following command :

set GDAL_DATA=D:\tests\test_qgis_import\dist\App\gdal
start App.exe

At this point, the application works correctly but is very slow, does anyone have an idea how to resolve this last issue ?

My python code :

import os
import os.path
import sys
from untitled import *
from mapTools import *
from PyQt4 import QtNetwork, QtXml, Qsci

class Main(Ui_MainWindow, QtGui.QMainWindow):

    def __init__(self, parent=None):...

    def EditPolygon (self):...

if __name__ == '__main__':

    QgsApplication.initQgis()
    app = QApplication(sys.argv)

    app.processEvents()
    frame = Main()
    frame.show()
    frame.raise_()

app.exec_()


texts: How to count pixels in a raster below a given value in PostGIS

texts:
The code below works in getting the basic statistics of a clipped raster. I would like to add a metric that is the percent of pixels that are below a given value. I am already calculating the total number of valid pixels, but I am missing the number of pixels below a given value. I have tried using ST_ValueCount, however this seems to only give the count of unique values in a band and not the count of all pixels below a given value. Does anyone know of a solution to my problem?

WITH
--  Select Features
    feat AS (SELECT DISTINCT fid AS building_id, wkb_geometry AS geom FROM inside_raster
        WHERE fid = 'osgb1000005284805'
    ),
--  Clip band 1 of raster image to the boundaries of buildings
--  then get statistics for these clipped regions
    b_stats AS
    (SELECT  building_id, (stats).*
        FROM (SELECT building_id, ST_SummaryStats(ST_Clip(a.rast,1,geom,-9999,true)) AS stats
            FROM height_raster a
                INNER JOIN feat
            ON ST_Intersects(feat.geom,rast)
        ) AS foo
    )
--  Summarise statistics
    SELECT building_id, count As pixel_count,
        min AS min_pixel_val,
        max AS max_pixel_val,
        mean AS avg_pixel_val,
        stddev AS pixel_stddev
    FROM b_stats
    WHERE count > 0;


texts: How can I create a fill beneath my 3D surface in ArcScene?

texts:
I came across this image on http://desktop.arcgis.com/en/arcmap/latest/extensions/3d-analyst/choosing-the-3d-display-environment.htm

I would like to create something similar that shows the below ground underneath my 3D TIN surface so that it looks like it's been sliced out.

However, I can't seem to find instructions for how to do this, what it's called or a tool that can generate it.

I've tried a few things such as the "extrude between" tool but I don't have a TIN to represent the bottom of the block to extrude between. I only have a TIN for the surface. I do have a polygon feature with the z value of 1m which is the exact extent of the TIN surface. I could use that as the base to extrude between but only get a 999 error when I try to generate TIN from it, probably because it is a flat plane rather than an irregular surface.

Are there any tools or work arounds for creating this?

I know it can be done in non-spatial software like InDesign if I export the surface and do some graphic design on it but my main skill in in GIS and it gets expensive and time consuming to work across multiple software packages.


texts:
I am building a map application, calling layers from an arcgis-server using the javascript api.

Maybe I am being stubborn, but I am trying to minimize my use of dojo throughout this process.  Much of my previous experience has been using other APIs like cartoDb, google maps and I prefer using jquery and pure javaScript.

I am going to build a dynamic list of layers and legend items by iterating through the layerInfos of the dynamicMapServiceLayer that I have added to my map.

I have found that I cannot itterate through the layerInfos until after the layer has been added to the map.  To detect this I am using the onLayerAdd event.  The only way I've been able to use the onLayerAdd event is to also use a dojo.connect

dojo.connect(map, "onLayerAdd",testFunction);

Do I have to use dojo.connect?  If so why?  It seems like it would make more sense to do something like

map.onLayerAdd(testFunction)

Thanks in advance....


texts:
Often I look at a map and wonder when it was created.

How do I figure out the age of an undated world map?

If I have a vector file of world political boundaries, but I am not sure of the date, then what resources can I use to determine the date?

Is there a map service on the web that allows me to enter a date and get a map of world boundaries for that date?

Assuming there is, what process would I use to determine how close the match is?

Are there any rule based systems that are specific to map dates?


texts: Is there a reverse geocoding service that lets me specify a lat/long bounding box?

texts:
I'd like to have all the street addresses within a given area, but I'm unable to determine if this is even possible. I can find many online services that given a specific latitude and longitude will return a single street address. The area I want this for in particular has approximately 500 residences, it's extremely tedious to attempt to compile the list by hand (and afterward, I may want another similar list later, I haven't decided).
The streets in question run parallel to north/south and east/west, so supplying the bounding box by either named streets or via latitude/longitude is simple enough.
Is there any easy way to obtain this data?


texts:
A reverse clip saves only the part of your spatial object that is outside
the bounds of another object, as opposed to a regular clip which saves the
parts that are inside the other object.

Performing reverse clip in ArcMap? shows how to do it in ArcMap.

How do I do this in R?

Reproducible example (on Linux machines):

system("wget 'https://github.com/Robinlovelace/Creating-maps-in-R/archive/master.zip' -P /tmp/")
unzip("/tmp/master.zip", exdir = "/tmp/master")
uk <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "ukbord")
lnd <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "LondonBoroughs")
plot(uk)
plot(lnd, add = T, col = "black")

What I want here to do is to save all of the UK except for London. Visually, I want the black shape in the resulting image to be a hole.


texts:
I'm trying to use Pyosmium to parse changesets from an osmChange file (to replicate the osmium changeset-filter [1] functionalities):

      class MyHandler(osmium.SimpleHandler):
      ...     def changeset(self,changeset):
      ...             print changeset

      handler = MyHandler()
      handler.apply_file('mychange.xml>')

but nothing happens. No errors, no outputs.

I've tried to setup the osmium.io.Reader manually, and setting the entity bit:

      reader = osmium.io.Reader('mychange.xml',osmium.osm.osm_entity_bits.CHANGESET)
      omsium.apply(reader,handler)

No output again.

Is it possible to parse osmChanges with Pyosmium? If yes, what's the right way to do it?

[1] http://osmcode.org/osmium-tool/manual.html#working-with-changesets


texts: iPad mapping application with custom maps

texts:
I'd like to load custom maps onto my iPad and use them offline with the GPS for navigation. These are maps I can prepare in any format and projection.

Is there an application that can take, for example, a GeoTIFF or group of GeoTIFFs and use them as the background with a GPS dot overlay?

Use case: I'm navigating about 100 km by snow-mobile across part of Antarctica. Normally we use standard GPS, but those screens are tiny, hard to see through goggles and hoods, and hard to interact with through gloves. I'd like to use an iPad (and a stylus to get around the gloved issue) because it has a nice big screen. I'll pack it with hand-warmers to keep it functional.

Data set: I have 0.5 m/pixel images of the terrain. At that resolution we can see the (visible) crevasses and the sensors (GPS, seismic) that we are trying to find, although they have moved a few 10s to 100s of m with the ice as it flows. Each image covers about 17 square km and is about 4 GB in size. I can generate these images as more smaller files, move up to 1 m/pixel if necessary, etc.


texts: Putting raster values in polygons ArcGIS

texts:
I'm trying to extract values from a raster file dataset (called land_use) to polygons (in a file called Hydro_buff) that I have created.
I read a few questions on that topic but some other issues come up.

I have tried "zonal statistics" but I get a raster output when I just need a feature output. So I tried using "zonal statistics as table" to then join the two tables together (join to Hydro_buff) but it doesn't work because for some reason I go from 74 034 polygons to just 31 426 after using zonal statistics as table so  lot of the polygons have Null value. All of my polygons have a UID that I created using the "create Unique ID" tool.

I checked the extent of both of my files and only 10 000 polygons don't overlap the land_use

Does anyone know what the problem is ?


texts: Can Polygon cross other Polygon on the surface?

texts:
I'm using Shapely (https://shapely.readthedocs.io/en/latest/manual.html)

I'm trying to get 'True' value for crosses method (https://shapely.readthedocs.io/en/stable/manual.html#object.crosses) with two polygons, but I can't find suitable data.

A have two polygons, that have an interior that intersects, The interior doesn't contain each other.   I have trying two cases. In the first case intersection was a polygon(later I learned that it was not suitable, the dimension of a polygon is 2). In the second case intersection was a boundary of the interior of two polygons(dimension 1). But the method crosses return False again. So I started to think that this check is not at all possible for two polygons on the surface


texts:
I am trying to install r.inund.fluv on Linux to produce inundation maps but I am getting the below error. Note that I have a Fortran compiler.

GRASS 6.4.4 (newLocation):~ > g.extension r.inund.fluv

Fetching <r.inund.fluv> from GRASS-Addons SVN (be patient)...
A    r.inund.fluv/r.inund.fluv_62
A    r.inund.fluv/correction_from_path.f90
A    r.inund.fluv/clean_inundation.f90
A    r.inund.fluv/r.inund.fluv
A    r.inund.fluv/find_main_channel.f90
A    r.inund.fluv/README.txt
A    r.inund.fluv/r.inund.fluv.html
A    r.inund.fluv/Makefile
A    r.inund.fluv/2d_path.f90
Checked out revision 65521.
Compiling <r.inund.fluv>...
Makefile:4: /usr/lib64/grass-6.4.4/include/Make/Script.make: No such file or directory
Makefile:5: /usr/lib64/grass-6.4.4/include/Make/Module.make: No such file or directory
make: *** No rule to make target `/usr/lib64/grass-6.4.4/include/Make/Module.make'.  Stop.
ERROR: Compilation failed, sorry. Please check above error messages.


texts: Difference between UTM CRS with 6 digits and 8 digits

texts:
In my area in Germany we use the CRS with the EPSG:25832 code. These are UTM 32N coordinates with a false easting of 500,000. Thus, it is a 6 digit coordinate system. I read this means that the middle meridian has the value 500,000m.
Sometimes I get data with the EPSG:5652 code which is also a UTM 32N System. But it has a false easting of 32,500,000. The false easting has 8 digits. I read that the 32 refers to the Zone 32.
What are the differences between both CRS?
In what situations is it recomended to use which version?


texts:
I am working with OpenLayers. I want to convert Shapefile into GeoJSON format and send it to browser to show with OpenLayers. But the problem is that when the size of my Shapefile is bigger than 60Mb, then it is not convenient to send whole GeoJSON data to browser.

Anyone can help me, how can i minimize the GeoJSON file size while my Shapefile size is still very large.

Are there any tools or techniques by which I can encode GeoJSON data to minimize the size and send it to client side as well as decode this data to show in map.


texts: Spatial Join by group in ArcGIS Desktop?

texts:
I have a point feature and a polygon feature and want to join the point to the polygon so that I can get the count for each polygon. Question is the points are grouped by a filed called type, it has like "A","B","C", etc.

Is there a way to join grouped points by the "type" field so that after joining, I'm going to have the polygon feature with fields like "CountA", "CountB", etc., and the it reflects the number of the points falls in that polygon respectively?

I don't see a query statement option in the spatial join function and currently my solution is to create different point features by "type" field, and make the spatial join respectively.

Is there an easier way to do that?

Polygon and points are stored in file geodatabase as below:

X:\DATA\DATABASE1.GDB\MAIN\POLYGON

X:\DATA\DATABASE2.GDB\POINTS1

Either in python window or stand-alone script will be fine.


texts: Maxent SDM projections coming up almost entirely blue

texts:
I have been generating SDM models and projections in Maxent. I have had some success with one set of variables (some Hydroshed variables mixed with World CLIM variables); however, I am having a bit of trouble with a model based solely on World CLIM variables.
I used a Pearson correlation matrix to select eight World CLIM variables, which model well with my species' occurrence data. The issue arises, when I ask Maxent to project on future World CLIM variables - the returned projection is almost entirely blue. There are a few specs of red and yellow (indicating prediction probabilities), but the vast majority of the projection is blue (indicating zero prediction probability). I tried rerunning the model with clamping turned off, but the problem persisted.
I am relatively novice with Maxent, so I was wondering if I am overlooking a setting that may help correct this issue. I am a bit confused, because when I ran the projection with the Hydroshed variables and some of the same World CLIM variables, the future projection looked just fine.
I am using ArcGIS Pro from Esri and the Maxent application (version 3.4.4) to make my models. I am taking rasters (ASCII files) of each variable, clipping them to my study extent, then feeding them through maxent with my occurrence data to generate the models.
I am attaching pictures of my current Maxent settings:

I am also attaching a picture of my functioning model (top map in the picture) and the problematic blue projection of future variables (bottom map in the picture):

As you can see in the future projection, there is some light wisps of color in the far North-East, so it appears the projection is at least being made without any major errors.
As a side note, I've played around with selectively using only certain regions of the occurrence data for modeling and noticed that certain selections will generate more normal looking future projections, while others will produce mostly blue maps, like the one in Picture 3.
Why are some of my future projections for only CLIM variables coming back almost completely blue, while projections made with CLIM and Hydroshed variables are coming back normal looking?


texts: I sent a crew out to collect data and take photos using Survey 123. The process worked great, but the crew used cell phones and took some photos with different orientations. Survey 123 doesn't reorganize and reorient the photos.  Does anyone know if you can edit the orientation of the attached photos without downloading the photos and uploading back to AGOL? I wanted to create a story map from the photos but don't see an easy way to do this.


texts:
I have an NDVI-Time series (24 tiffs from one year). What I would like to calculate is one single Raster with the Maximum NDVI-value for each pixel from all 24 images. I thought about using gdalcalc but I cannot find the right syntax to calculate the maximum of a my inputfiles:

Example for 2 Layers:

shell("E:/Programme/OSGeo4W.bat; gdal_calc -A E:/IVFL/Somalia/NDVI/201401_001_NDVI_NAfill.tif -B E:/IVFL/Somalia/NDVI/201409_121_NDVI_NAfill.tif --outfile E:/IVFL/Somalia/NDVI_maximum.tif --calc=(max(A,B))")

calc=max(A,B) always gives me an error.

And in the end I even have 24 input files and not only 2.
Does anyone know the syntax (numpy) to calculate the maximum?


texts: The total area of all land cover data exceeds the corresponding polygon area

texts:
I am using the national land cover raster data. I have an MSA boundary shapefile.  I want to calculate the area of each land cover within each MSA area. So I used the Tabulate Area tool in ArcGIS Pro to calculate it. But I found for each MSA, the total area of all land cover exceeds the MSA area.
Can someone help me with it?
The resolution of the MSA shapefile is 1:500,000.
the resolution of land cover data is 30 meters.
I have set the cell size to be consistent with the land cover data. But still, I have this problem.


texts: Where to get the modern Constellations as a vector File?

texts:
Does anyone know where I can get a vector version (shapefile/TAB) of the modern star constellations for the entire sky?

Wikipedia has a SVG - https://en.wikipedia.org/wiki/File:Constellations_ecliptic_equirectangular_plot.svg - but I'd prefer not to have to try and convert it to a "proper" format.

Basically I want the below, but as a vector; ideally with the attribute data of constellation name for each polygon.

Thanks.


texts: How to crop a shapefile using a raster with Geotools?

texts:
I'm working on a project for my school for which I need to use Geotools.
I need to be able to display a raster file and a shapefile, but the shapefile is way bigger than the raster and I just want to display the part in common with the raster.

Would someone know how to "crop" my shapefile and only keep the vectors in a specific area ? With Geotools of course.
I'm using a Tiff with a tfw file as a raster.


texts: MYSQL 5.6 Cluster / chain points by distance to any other point without millions of calculations?

texts:
I have a table of approx 10,000 - 20,000 dynamic points across the globe stored as lon/lat points. I want to cluster any point within (roughly) 3km of any other point into one point (and return the averaged center point of each cluster if possible). Unlike some other solutions I have seen, I do not just want points within 3km of the central cluster point, it is fine to have a point within each cluster many kilometres away from the central cluster point as long as there is a chain of points within 3km of each other joining the points together into a single chain/cluster.

The only way I can think of doing this is looping over every point, finding any other points within 3km of that point, and then recursively each of those points etc until i find all in that perticular chain/cluster, and then repeat for all the remaining un-clustered points - this would result in millions of calculations. Is there an easier way of doing this? Preferably with MYSQL 5.6 geospatial abilities, but im open to any viable solution.


texts: Setting Buffer from GeoprocessingTools with a feature

texts:
I'm trying to set up a Buffer from Geoprocessing Tools and i would like to put a Field from the layer as Distance Souce.

Here is my code :

import processing
layer =  QgsVectorLayer(path_to_routes, 'Routes', "ogr")
QgsProject.instance().addMapLayer(layer)
processing.run('native:buffer', {"INPUT": layer, "DISTANCE":50, "OUTPUT": memory:buffer})

I don't find a solution to put a field on "DISTANCE" parameter.

Here is a picture of what i want to set up.


texts: Display arrows on each segment of a multines feature, using GeoServer and SLD

texts:
I'm using the SLD code below to display arrows along lines with GeoServer 2.10:

<Rule>
  <Name>Lines</Name>
  <Title>Lines</Title>
  <LineSymbolizer>
    <Stroke>
      <CssParameter name="stroke">#000000</CssParameter>
      <CssParameter name="stroke-width">1</CssParameter>
    </Stroke>
  </LineSymbolizer>
  <PointSymbolizer>
    <Geometry>
      <ogc:Function name="centroid">
        <ogc:PropertyName>the_geom</ogc:PropertyName>
      </ogc:Function>
    </Geometry>
    <Graphic>
      <Mark>
        <WellKnownName>shape://oarrow</WellKnownName>
        <Fill>
          <CssParameter name="fill-opacity">0</CssParameter>
        </Fill>
        <Stroke>
          <CssParameter name="stroke-width">1</CssParameter>
        </Stroke>
      </Mark>
      <Size>20</Size>
      <Rotation>
        <ogc:Function name="endAngle">
          <ogc:PropertyName>the_geom</ogc:PropertyName>
        </ogc:Function>
      </Rotation>
    </Graphic>
  </PointSymbolizer>
</Rule>

This is doing the job properly when features are single lines (straight segments). In the case a feature is a multiline including segments going in various directions, the arrow is displaying wrongly, taking into account the segments all together (on the image below, the arrow is cercled in blue).

The solution would be to explode the multilines in the database table on which the GeoServer layer relies, but I need to keep my features as multilines. Any idea on how to achieve that?


texts: Extract Polygon Coordinates from Shapefile data

texts: I have some building data (polygons) as a shapefile in QGIS 2.18. I'm trying to add a location coordinates (the polygons rough center) and an array of the polygons coordinates (its points). I've been able to do the first one using vector --> analysis tools --> mean coordinates which generates a new layer. I can't work out how to get the polygons coordinates though? The end goal is to extract this all to a spreadsheet if that makes a difference.


texts: Looking for a Tool, that allows me to "Union" many overlapping layers, while prioritizing some to overwrite the others in ArcGIS

texts:
I am currently working on a project where I analyze the contents of a 500m buffer around several sample points.
Foundation for those informations are multiple shape file layers, which, alltogether, cover the area without having blank spaces. Unfortunately, this also results in overlapping as seen here.

Due to different dates of creation of the data, I want to prioritize recent data and not older ones. So during a Union of my data, I would like to know if there is a way for hand-selected, "newer" data to overwrite older data, thus resulting in a more compromised data table which would then enable me to pie chart the contents, according to the percentage of area it has in my 500m buffer Clip.
I would also accept how to get to my goal differently than by messing with Union and Table content managing.


texts: Make a series of maps showing different rasters

texts:
I have a set of rasters and I would like to create a series of maps, in which each page/map shows a different raster. The spatial extent is the same among all rasters. How could you achieve this in ArcGIS Pro?

I have looked into using the map series option, however it is focussed on the same items on the map but different spatial extents. I thought about making the map series index like a dummy: many pages which would focus on an identical area of interest. I know from ArcMap that it's possible to do page queries to only display data if a property matched with the data driven page. However, this only works for not-raster-features. I also stumbled upon Raster Catalogs which seemed to have potential but that is no longer supported in ArcGIS Pro.

Does anyone have an idea how I could achieve this in ArcGIS Pro?


texts: How do I set up a link to a WMSLayer in JOSM? URL works in QGIS but not JOSM

texts:
I have a publicly available WMS link that I'm trying to use as a layer in JOSM. It comes from http://dggs.alaska.gov/

http://maps.dggs.alaska.gov/arcgis/services/apps/Lidar_Hillshade/MapServer/WMSServer

When I load it in QGIS it connects to the server and there is a list of layers that I can add. When I try to load it in JOSM, there are no options of choosing layers, just inputting the link. I use 'add rectified layer' and enter the URL but no tiles are displayed.


texts:
I am wondering if there is a way to show an area of a polygon as you draw it. I trying to create a perfect 1 acre polygon, but the polygon could have 1-20 sides. I don't really want to have to do $area in the attribute table every time to find my area.
Any suggestions?


texts: R package gapfill: How to convert R raster stack to 4-dimensional array and then back to raster stack?

texts:
I'm currently analyzing MODIS EVI time series (MOD13Q1) for a small region between 2000 and 2015. Now, I would like to fill missing values using the new R package gapfill.

What I have done so far:

Downloaded and pre-processed data using R MODIS package. This gave me 377 GeoTiff-files that are named MOD13Q1.A2000081.250m_16_days_EVI.tif and so on.
Used the QA layers from MOD13Q1 to remove pixels with clouds.

As a result, I now have 377 GeoTiffs with NA values for clouds:

> stack(evi_without_clouds_file_paths)
class       : RasterStack
dimensions  : 121, 122, 14762, 377  (nrow, ncol, ncell, nlayers)
resolution  : 257.2911, 257.3028  (x, y)
extent      : 448840.3, 480229.9, 7766995, 7798129  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=38 +south +datum=WGS84 +units=m +no_defs     +ellps=WGS84 +towgs84=0,0,0
names       : MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, ...
min values  :                 -1195,                  -846,                  -220,                  -185,                  -321,                  -344,                  -128,                  -159,                  -309,                    99,                  -281,                  -268,                  -200,                   218,                  -290, ...
max values  :                  9205,                  9902,                  8466,                  7933,                  7696,                  6700,                  7019,                  6825,                  5578,                  5163,                  4934,                  5028,                  4561,                  5691,                  4262, ...

However, the input data for gapfill has to be in the following format: "Numeric array with four dimensions. [...]. the data should have the dimensions: x coordinate, y coordinate, seasonal index (e.g., day of the year), and year."

I'm familiar with the raster package, but I have no experience whatsoever working with multidimensional arrays. As a result, I spent the whole day trying to convert my data, but had no success so far. So my questions would be:

How to convert raster stacks to 4-dimensional named arrays?
How to convert 4-dimensional named arrays back to either a brick with 377 layers or 377 single GeoTiff files after having gaps filled with gapfill?

I apologize for not being able to provide sample data. I tried to simulate similar data, but didn't manage to do so.

I know that it is difficult to answer the question without sample data, but I would also be very happy about some pointers to the right direction.

Here is the closest I have gotten so far:

require(tidyr)
evi.stack <- stack(evi_without_clouds_file_paths)
x <- as.data.frame(evi.stack, xy = TRUE)
x <- gather(x, key = "scene", "value", -x, -y)
x$scene <- extract.id(x$scene)
x$day.of.year <- as.numeric(str_sub(x$scene, 6, 8))
x$year <- str_sub(x$scene, 2, 5)

evi.dim.long <- unique(x$x)
evi.dim.lat <- unique(x$y)
evi.dim.day.of.year <- unique(x$day.of.year)
evi.dim.year <- unique(x$year)

x.vec <- x$value
dim(x.vec) <- c(length(evi.dim.long),
                length(evi.dim.lat),
                length(evi.dim.day.of.year),
                length(evi.dim.year))
dimnames(x.vec) <- list(evi.dim.long,
                        evi.dim.lat,
                        evi.dim.day.of.year,
                        evi.dim.year)

If I only import data from a single year, this seems to kind of work (although the resulting image was "upside down"). I was able to fill gaps with Gapfill, but then had no idea how to convert data back to a format that could be further processed with the raster package.
If I try to import data from more than one year, the dim command fails.


texts: I've been creating tools for ArcMap using Delphi through COM even though it isn't officially supported and I'd like to know if that will still be possible with ArcGIS Pro?


texts: In QGIS for Type Case in Labels, Rule Based Rule, the Capitalize first letter option does not work

texts:
Here are my labels untouched

Now when i change to all lowercase

You see that it changes

However for type case i want to capitalize first letter

It shows that it made it everything capitalized

The default lettering of the values are all capital.

I tried the expression below. "Street" is the name of the field. It even acknowledged the output preview and changed it how i would have liked to ideally.

title("Street")


texts: Creating a openlayers enabled web project in eclipse

texts:
I'm brand new to  creating gis map applications and I am trying to get a web project that uses

function init(){
            map = new OpenLayers.Map( 'map' );
            layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
                    "http://vmap0.tiles.osgeo.org/wms/vmap0",
                    {layers: 'basic'} );
            map.addLayer(layer);
            map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);
            var featurecollection = {
              "type": "FeatureCollection",
              "features": [
                {"geometry": {
                    "type": "GeometryCollection",
                    "geometries": [
                        {
                            "type": "LineString",
                            "coordinates":
                                [[11.0878902207, 45.1602390564],
                                [15.01953125, 48.1298828125]]
                        },
                        {
                            "type": "Polygon",
                            "coordinates":
                                [[[11.0878902207, 45.1602390564],
                                  [14.931640625, 40.9228515625],
                                  [0.8251953125, 41.0986328125],
                                  [7.63671875, 48.96484375],
                                  [11.0878902207, 45.1602390564]]]
                        },
                        {
                            "type":"Point",
                            "coordinates":[15.87646484375, 44.1748046875]
                        }
                    ]
                },
                "type": "Feature",
                "properties": {}}
              ]
           };
           var geojson_format = new OpenLayers.Format.GeoJSON();
           var vector_layer = new OpenLayers.Layer.Vector();
           map.addLayer(vector_layer);
           vector_layer.addFeatures(geojson_format.read(featurecollection));

        }

I've included the OpenLayers.js with the img and theme folders under the WEB-IMF.
When I run it only the text shows, not the map.
Is there something I'm missing.


texts:
I got this error when running osmosis:

root@osm:/home/osm/updates# osmosis --rri workingDirectory=. --wxc ./update.osm.gz
Apr 13, 2017 10:03:34 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Osmosis Version 0.45
Apr 13, 2017 10:03:35 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Preparing pipeline.
Apr 13, 2017 10:03:35 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Launching pipeline execution.
Apr 13, 2017 10:03:35 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Pipeline executing, waiting for completion.
Apr 13, 2017 10:03:37 AM org.openstreetmap.osmosis.core.pipeline.common.ActiveTaskManager waitForCompletion
SEVERE: Thread for task 1-rri failed
org.openstreetmap.osmosis.core.OsmosisRuntimeException: Unable to read the changeset file 000/001/672.osc.gz from the server.
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.downloadReplicationFile(BaseReplicationDownloader.java:135)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.download(BaseReplicationDownloader.java:264)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.runImpl(BaseReplicationDownloader.java:304)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.run(BaseReplicationDownloader.java:383)
        at java.lang.Thread.run(Thread.java:745)
Caused by: java.io.IOException: No such file or directory
        at java.io.UnixFileSystem.createFileExclusively(Native Method)
        at java.io.File.createTempFile(File.java:2001)
        at java.io.File.createTempFile(File.java:2047)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.downloadReplicationFile(BaseReplicationDownloader.java:113)
        ... 4 more

Apr 13, 2017 10:03:37 AM org.openstreetmap.osmosis.core.Osmosis main
SEVERE: Execution aborted.
org.openstreetmap.osmosis.core.OsmosisRuntimeException: One or more tasks failed.
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.waitForCompletion(Pipeline.java:146)
        at org.openstreetmap.osmosis.core.Osmosis.run(Osmosis.java:92)
        at org.openstreetmap.osmosis.core.Osmosis.main(Osmosis.java:37)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:57)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:606)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launchStandard(Launcher.java:330)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:238)
        at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:415)
        at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:356)
        at org.codehaus.classworlds.Launcher.main(Launcher.java:47)

My configuration.txt :

# The URL of the directory containing change files.
baseUrl=http://ftp5.gwdg.de/pub/misc/openstreetmap/planet.openstreetmap.org/replication/day/

# Defines the maximum time interval in seconds to download in a single invocation.
# Setting to 0 disables this feature.
maxInterval = 86400

I CAN get the same file using wget:

http://ftp5.gwdg.de/pub/misc/openstreetmap/planet.openstreetmap.org/replication/day/000/001/672.osc.gz

Behind a proxy, but this not solve:

export JAVACMD_OPTIONS="-Djava.io.tmpdir=/tmp -Dhttp.proxyHost=xxx.xxx.200.10 -Dhttp.proxyPort=999 -Dhttps.proxyHost=xxx.xxx.200.10 -Dhttps.proxyPort=999 -DsocksProxyHost=xxx.xxx.200.10 -DsocksProxyPort=999"

Some clues?


texts: Calculate volume from raster (GRASS) with depth per cell

texts:
I'm trying to calculate the volumes of depression areas filled with water. I have the depth for each cell in the raster, and I'm trying to get the total volume of cells adjacent to each other. Which I then, in step 2, need to remove all the bodies of water lower than a given threshold and then convert the file to vector/shp format.

I've trying using the r.clump and then r.volume function in qgis, but I keep getting this error:

ERROR: r.volume: Sorry, is not a valid parameter
ERROR: Required parameter not set:
(Name of input raster map representing data that will be summed within clumps)

C:\PROGRA~1\QGIS2~1.18\bin>v.out.ogr -s -e input=centroids0814b2718bf14943909f424cebd714fb type=auto output="C:\Users\JGJ\AppData\Local\Temp\processing503b826c304e474ebc5a9862a9dd29e1\b97bc604e110497ab8dde2cc3da478e7" format=ESRI_Shapefile output_layer=centroids
ERROR: Vector map not found

C:\PROGRA~1\QGIS2~1.18\bin>exit
Converting outputs
Loading resulting layers

The following layers were not correctly generated.
Centroids
You can check the log messages to find more information about the execution of the algorithm

I could not get the above to work, but managed a workaround.

I have converted raster to vector, and calculated each polygons volume. I have a lot of polygons I want to dissolve into a lot of larger polygons, which can be done with Dissolve. However each of the small polygon has a volume too it, that needs to be summed up. If I do a dissolve it just adds the total volume of all the polygons in the layer, and not those that get dissolved. Added a concept picture.


texts:
I am trying to update the FACILITYID field with a String "SWW" and then follow the "SWW" with the OBJECTID field.

So I will end up with SSW1", "SSW2" all the way to the end in the FACILITYID field.

I know for python it is simply: "SWW" + str(!OBJECTID!)

How would you do this using VB in field calculator?


texts: What operations do the GPS recievers perform after capturing signal from Satellite?

texts:
What are the operations on the signals recieved by a GPS reciever ?
How does it process them towards achieving a position fix ?

I need explanation with the decoding and correlating part of the signal. How does it determine the position of a satellite when the signal it is recieving from many satellite is mixed up with the signals from other satellites in the wireless medium .

I'd be glad if you point me to a good article/link that explains this well.
Thanks


texts: Using filters to filter a variable value with GeoTools

texts:
I have a Java Swing application which loads a line shapefile and draws a map. I'm using a filter and a rule to filter out the lines in the shapefile where the attribute "traffic" is greater than 500. Each line has a different "traffic" value and this works perfectly ("traffic" means the maximum allowable traffic volume of the road). I can make a separate style for this filtered lines.

    Rule rule = sb.createRule(textSymbolizer);
    rule.setFilter(sb.getFilterFactory().greater(
            sb.getFilterFactory().property("traffic"),
            sb.getFilterFactory().literal(500.0)));
    FeatureTypeStyle fts = sb.createFeatureTypeStyle("Feature",rule);

Now I want to do the following logic

if("traffic" property value > current_traffic_value){
    // style this line uniquely
}

The "current_traffic_value" is a integer value stored in a database for each line. So I want to check each line against a variable value. Can I achieve this using SLD methods? What would be the optimum method to achieve this?


texts: Exporting from Google Earth Engine into Google Earth Pro

texts:
I've found a dataset (Link Here) that I would like to export into Google Earth Pro so that I can analyse it alongside other data. Is this possible for me to do?
I'm sorry if the question has a simple answer but I have next to no experience with coding - I just need the dataset into a file that could be imported into Google Earth Pro. I can't seem to find any simple export button and a google search seems to suggest there's some way of coding an export? But this is totally beyond me at the moment.


texts: Find and replace text in metadata of feature classes

texts:
I have 36 feature classes in individual feature datasets.  I need to update a portion of the metadata for all the feature classes.

I want to find a string of text in the distribution portion of the metadata and replace it with a new string of text.

What would be the best way to do so? Do I have to manually update all 36 feature classes?

Can you write a python script to do that?


texts: Spatialite table column is function in itself

texts:
GOAL: I'm trying to filter all the 1's in a spatial boolean operation in Qspatialite, the spatialite plugin for QGIS.

ATTEMPTED SOLUTION: My solution is to perform the boolean query first, export the resulting table AS a table second, then query THAT table for all the 1's, third.

PROB W/ SOLUTION: Problem is, column to be filtered has a header that is in itself a function (i.e., "within(geom 1, geom 2)"). Thus, the header cannot be called in the command dialog as a simple string, because it inevitably becomes a command/function.

This is what I ran:

SELECT Towns.name, Schools.name, within(Schools.geom, Towns.geom) from Towns, Schools

So, within(Schools.geom,Towns.geom) becomes the table header, which I need to filter for 1s.

Any ideas/suggestions?


texts:
I was wondering if anyone knows of geospatial tools that work in SketchUp?

I know if GeoSketch, which allows a user to import ortho imagery, DEMs, shapefiles and GPS enabled photographs into SketchUp for 3D model creation.  I am looking for other geospatial tools that are designed for SketchUp, but I am having issues finding any.

I am interested in geospatial tools that work within SketchUp.


texts: Surface area calculation in Google Earth Engine?

texts:
Using this code, I've extracted water extent in a Landsat median image. But how can I compute water area extent?

code link: https://code.earthengine.google.com/8ad84d9ca59b3cb71b62f83119621bdc

var urmia = table;

Map.centerObject(urmia);
Map.addLayer(urmia);

var landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterBounds(urmia)
.filterDate('2018-01-01','2019-01-01')
.filter(ee.Filter.lessThan('CLOUD_COVER',5))
.median()
.clip(urmia)
.select('B[4-5]')
.normalizedDifference(['B4','B5'])
.gt(0.2)
.rename('time2');

Map.addLayer(landsat);


texts: Multiple Linear Regression Analysis in QGIS (with SAGA Toolbox)

texts:
I'm attempting to run a multiple linear regression analysis using QGIS, through the SAGA toolbox.

However, the QGIS UI does not allow me to choose predictors -- it seems entirely absent from the UI. Is this a bug, or am I mis-understanding something?

I am running QGIS 2.18.15, with SAGA 2.3.1 under Ubuntu 17.10.


texts: QGIS creating batch folder directory with incremented names by using Python console

texts:
I have the directory, where I have to store all the items separately in the list of folders. This list includes exactly the same folders only with incrementing number, as you can see below:

Because I am going to have more than 200 folders, I decided to create all of them by using Python console in QGIS.
My approach so far looks as you can see below:
 import os
 path='C:\\Users\\mk\\Desktop\\AX\\Pre survey\\Tyton'
 os.chdir(path)
 for i in range (1,208):
 Newfolder='Strand-' + str(i)
 os.makedirs(Newfolder)

As per the tutorial below:
https://www.youtube.com/watch?v=4nr0eW3w1u4&ab_channel=PyBTV
Since I managed to create one folder, there is a problem with the batch of these folders.
First of all, I cannot put something like "leading 0", as I tried:
    for i in range (001,208)

I cannot do this, but I need the folder numeration as shown in the image.
Another issue is, that if I use the whole code presented, Python throws an error:
IndentationError: expected an indented block
debugging the following line:
    Newfolder='Strand-' + str(i)

The examples how to solve this error are here:
https://stackoverflow.com/questions/4446366/why-am-i-getting-indentationerror-expected-an-indented-block
https://stackoverflow.com/questions/492387/indentationerror-unindent-does-not-match-any-outer-indentation-level
but i couldn't apply them to my situation.
How can I do the batch directory creation with Python, preferably by QGIS python console?


texts: Displaying Live Incoming Smartphone Images with Ushahidi on Map?

texts:
I have not yet used Ushahidi and I could not find the answer on their website to my question.

I wish to do the following:

setup an Ushahidi server/map,
go outside to make photos with my Android/iPhone smartphone
upload these photos to Ushahidi which should be saved in a database
display the photos on the map according to the coordinates. Also few lines of comments would be nice.

When I now visit my site and open the map, I would like to see some dots (actual position of an entry/uploaded image). By double clicking on the dot, the original uploaded image+comments should be displayed(pop up).

Is this possible with Ushahidi?


texts: Openlayers 2 How to know the map is after panning/zooming

texts:
Firstly i found a logic to customize or extends Openlayers 2 Control Navigation. I dont want to custom the openlayers default, because I'm currently using 2 different map in one page. especially to custom after panning or zooming.

The purpose is I Whenever user after panning/zooming/wheel-zoom or the position/zoom of first map is change, it can get latests 4 LatLon of corner side from frame viewed first map then call a function to process it.

edit:
already found a way using event register


texts: Historical OpenStreetMap data and dummy lines

texts:
I'm trying to export some historical highway data from OSM  with the following commands:

osmium time-filter clipped_greece_history-190715.pbf --overwrite  "${YEAR}-12-31T00:00:00Z" -o "osm${YEAR}.pbf"

# extract specific categories
osmosis \
    --read-pbf file="osm${YEAR}.pbf" \
    --tf accept-ways highway=motorway,trunk,primary,secondary,tertiary,unclassified,residential \
    --used-node \
    --write-pbf "osm${YEAR}_highways.pbf"

But I get some dummy(?) lines:

Am I missing something?
I cannot find an easy way to exclude those lines. Some of them are not straight lines.


texts: Interact with Selection and Entire Dataset using arcpy

texts:
ArcToolbox geoprocessing tools (and therefore arcpy) tend to run either on an entire dataset (if there is no selection) or on a selection set (if it exists).  How can I set up an arcpy script that will have some tools run on the entire dataset while others run on just the selection set?

Scenario - I have a script that does some geometric network traces based on a selection of water mains.  When this is run the script modifies the Enabled flag on a number of network points to modify the trace results.  I have found that I also need to modify some of the network lines, however when I run the tool that updates this flag it only runs on selected lines not on the entire dataset.

How can I make this portion of my script update lines anywhere in my network (whether selected or not) without changing the selection state of the currently selected records (needed for later in the script)?

I thought maybe I could do this with Make Feature Layer however this also uses just the currently selected features.

The only two ways I've come up with so far are:

to store the selection from the FIDset, clear the selection, make the modification, then re-select
Duplicate the layer in ArcMap and have the script reference both

Neither of these are optimal, and I think things would work better if I could just specify somehow (env setting maybe?) for arcpy to honour or ignore the selection set.


texts:
I am trying to access the WMS layers from the UK Coal Authority website. I can connect to their servers and see the raster files available but when I download them, I can only view the rasters from the first and last services:

Coal mining reporting areas - https://map.bgs.ac.uk/arcgis/services/CoalAuthority/coalauthority_coal_mining_reporting_areas/MapServer/WMSServer - the areas in which the Coal Authority recommend a mining report.
INSPIRE Compliant Data - https://map.bgs.ac.uk/arcgis/services/CoalAuthority/coalauthority_inspire/MapServer/WMSServer - dataset comprising Coal Authority monitoring points and licence areas.

The others do not work for me. These were my steps:

Create new WMS connection and only add a Name and copy URL location from website:

Click Connect and select a raster layer:

Layer is added along with a legend but nothing is shown on the map canvas. I also checked the properties of the layer and noticed that there were no rows/columns of data:

Is there something wrong with my procedure or with the rasters?

Doubt it's relevant but I am using QGIS 2.16.0-Nødebo.


texts:
L.Control.fileLayerLoad({
    fitBounds: true,
    layerOptions: {
        style: style,
        pointToLayer:function (data, latlng) {debugger;
            return L.circleMarker(
            latlng,
            { style: style }
            );
        },
        onEachFeature: function (feature, layer) {
            console.log("FEATURE",feature,"LAYER",layer);
            geojsonLayer = layer //
        }
    }

}).addTo(map);

I am trying to collect geojson data and use it. this function gives correct results for line and polygon. But when I use point data it gets stuck.


texts: Concatenating two fields (type double) into single field?

texts: In ArcGIS 10.3 using the field calculator, I'm trying to concatenate a census tract field (double) with a census block field (double) into a single field that has the two values separated by a decimal point (ex: 501.3001 ) and then generate a summary table that has counts for how many times each census block occurs. I'm having trouble combining these two fields with a decimal point in between them.  Should the new field be text? double? floating point?