texts: Using Mosaic to New Raster in ModelBuilder?

texts:
I have a large number of airphotos giving coverage of an area going back decades. - someone may have a smoother means to accomplish my goal here. I'm weak in Python sadly. I'm attempting to shrink/trim/eliminate a given number of pixels along the edges of the .TIF airphotos that I have.

The photos are quite old and scanned from negatives, and have black borders from the scanning of the film negatives, so I can't set nodata=0,0,0 as the border isn't a perfect 0,0,0 (in some cases water features are closer to 0,0,0 that then border is)

There's too many images to sensibly do this individually. I'll walk through my thoughts/plan below. My end goal is to use Build Footprints 'Shrink Distance' Feature:

Use recursive iteration to go through every .geotiff airphoto within the "airphotos" folder.
Create a file geodatabase where the name of the .gdb matches the raster
Create a Mosaic Dataset within that File Geodatabase [Mosaicmb]. The tools run as hoped up to this point, i've been able to create all the geodatabases.
Here the problems start: Add the raster that is being iterated to the newly created empty mosaic dataset. Fails with the error listed below
Build Footprints, and use the shrink feature to eliminate my borders, by cutting say 50 pixels off all sides of each photo.

ArcGIS reports that mosaicmb doesn't exist - which it obviously doesn't, and it won't exist until the tool runs.

If you know a slick way to trim border pixels, that will help too.

Error 000732 - Mosaic Dataset: Dataset Q:\.......\test5\%name%.gdb\mosaicmb does not exist or is not supported

I'm still experimenting and trying what works, but at this point I'm shooting in the dark and praying. I'll update if I have a breakthrough.


texts: GeoDjango / PostGIS - Cut geometry by bounding box

texts:
I am looking for a solution to retrieve the geometry(ies) from a polygon/multipolygon within a bounding box in a GeoDjango application.

The application I am developing has a requirement to serve out tiles of geographical data, and many of our geometries are very large and span across multiple tiles. To reduce the payload of each tile I would like to splice(cut) a geometry shape(s) from a bounding box.

I am storing Planet.OSM data in a PostGIS database.

Any suggestions would be greatly appreciated!

Thanks


texts: How can I fix a geoTiff with an incorrect CRS?

texts:
Trying to load and plot a raster (it's a DEM) in R, but having problems.

Using R I read in a raster file:

library(raster)

dem <- raster("DEM_Victoria.tif")  # >70 MB
print(dem)

class       : RasterLayer
band        : 1  (of  3  bands)
dimensions  : 5444, 6007, 32702108  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent      : 176615.6, 777315.6, 5650397, 6194797  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +no_defs
data source : DEM_Victoria.tif
names       : DEM_Victoria
values      : 0, 255  (min, max)

My confusion here is the extent vs the CRS. The extent shows coordinates from a UTM Zone 55S system, while the CRS reports it as longlat.

Similarly, gdalinfo reports:

Driver: GTiff/GeoTIFF
Files:  DEM_Victoria.tif
Size is 6007, 5444
Coordinate System is:
GEOGCS["GDA94",
    DATUM["unknown",
        SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (176615.587999997427687,6194797.138000000268221)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  176615.588, 6194797.138) (Invalid angle,Invalid angle)
Lower Left  (  176615.588, 5650397.138) (Invalid angle,Invalid angle)
Upper Right (  777315.588, 6194797.138) (Invalid angle,Invalid angle)
Lower Right (  777315.588, 5650397.138) (Invalid angle,Invalid angle)
Center      (  476965.588, 5922597.138) (Invalid angle,Invalid angle)
Band 1 Block=6007x1 Type=Byte, ColorInterp=Red
Band 2 Block=6007x1 Type=Byte, ColorInterp=Green
Band 3 Block=6007x1 Type=Byte, ColorInterp=Blue

Is this real? A problem? (well it is to me, but perhaps because I don't understand something) Fixable?

I tried 'Assign projection' in QGIS, which uses this command:

gdalwarp -t_srs EPSG:28355 -srcnodata -32768  DEM_Victoria.tif DEM_Victoria.tif.tmp

but it failed with:

  ERROR 1: Too many points (12100 out of 12100) failed to transform,
  unable to compute output bounds.


texts: Python console not loading with QGIS1.7-Ubuntu 10.04.2 LTS lucid,

texts:
Any tips on how to get the Python console working on Q-GIS (Quantum GIS)? (Its not working from the software's menu.)

System:

Python 2.6.5
QGIS1.7 OS
Linux: Ubuntu 10.04.2 LTS lucid

Error message from python running in the Terminal window:

Python 2.6.5 (r265:79063, Apr 16 2010, 13:09:56)
    [GCC 4.4.3] on linux2
    Type "help", "copyright", "credits" or "license" for more information.
    >>> from qgis.core import *
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
    RuntimeError: the sip module implements API v7.0 to v7.1 but the qgis.core module requires API v6.0
    >>>


texts: Why won't data driven map center and crop each township?

texts:
I need to create maps in ArcGIS Desktop 10.7.1 for many townships which are:

centered on each township
extent cropped at the border
fixed scaled at 1:24,000.

Fig. 1 shows what I would like. However, as I look through various townships they do not all line up nicely in the frame. Fig. 2 is too far south and the extent of Fig. 3 is too large (red marks indicate the desired extent corners). Fig. 4 and 5 show the settings that I have for the data driven pages.


texts:
I wanted to re-project a QGIS project, but there is no such option. When searching for an answer I see ogr2ogr suggested for this task, by re-projecting each vector layer (shapefile in this example). However I don't get expected result:

C:\temp>ogrinfo -al -so misc.shp
INFO: Open of `misc.shp'
      using driver `ESRI Shapefile' successful.

Layer name: misc
Geometry: Line String
Feature Count: 10
Extent: (21.267388, 42.015857) - (21.270225, 42.017470)
Layer SRS WKT:
PROJCS["MGI_Balkans_zone_7_deprecated",
    GEOGCS["GCS_MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",21],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",7500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
id: Integer (10.0)

C:\temp>ogr2ogr -t_srs EPSG:4326 misc_re.shp misc.shp

C:\temp>ogrinfo -al -so misc_re.shp
INFO: Open of `misc_re.shp'
      using driver `ESRI Shapefile' successful.

Layer name: misc_re
Geometry: Line String
Feature Count: 10
Extent: (-34.004490, 0.000230) - (-34.004490, 0.000230)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
id: Integer (10.0)

So new file has changed geographic coordinate system, but it's not projected.
How to do that?


texts: Is there any plugin that plots a profile from point survey data in QGIS 3

texts:
I used to work with the plugin ProfileFromPoints, but it is not avaible 8or updated) in Qgis 3.

I wonder if there is another one which does the same or if there is an easy (foolproof) way to use the same plugin in Qgis 3.


texts: Display and style a single US state in Mapbox/OpenStreet?

texts:
Ok, be gentle. I'm self-taught and don't know the 100% technically correct way to phrase my question.. :)

My goal is to create a Mapbox style that shows only the state of Virginia, and any standard features like roads, place names, water that I want to add. But only show it for Virginia.

Within Mapbox I have created a style that shows the state of Virginia (and only the state of Virginia) by nature of importing a JSON file of all the counties in the US, then filtering that layer to show only counties in Virginia. But adding in a layer like place names adds them across the entire US and they show up on the map I want restricted to Virginia. I don't see a state filter to help with some of this additional layer data.

It was suggested to use Overpass Turbo to query only what I need, but I don't see a way to limit a query to just Virginia.

Any suggestions (written in novice terms) to get my map to only show a Virginia map?


texts:
I would like mask several images in bash shell, gdalwarp works for me as a simple single line command:

gdalwarp -cutline mask.shp -cwhere 'kepnev="a.tif"' -crop_to_cutline a.tif a_masked.tif

..but it dosen't work with for loop:

for kep in *.tif;do
    gdalwarp -cutline mask.shp -cwhere 'kepnev="$kep"' -crop_to_cutline $kep ${kep/.tif}_masked.tif;done

Please help. What am I doing wrong?


texts: ReClassify Vector Data Between Values using ArcGIS Desktop?

texts:
I would like to be able to take all values between 1 and 10 and have a field define them as 1, and all values between 10 and 20 and set them to 2, how do I do that with vector data?

I read how to set all building to 20, and grass to 40, but I'm pretty sure that to do that for a range of values, like 1-10, I would have to manually enter each number, and decimal in that range. Am I wrong? and if not, how would I accomplish this?


texts: A Fast Way to Color Polygon With 200 Row With Six Selected Color in QGIS 3.4.2

texts:
I have a  layer of the world that contains 200 rows.

I am required to color the map with six selected colors, the rainbow colors.

I figure that I can save selected countries and save them separately.

Once they saved I can upload them and make the color ramp I needed.

However, I was wondered is there a faster way to do this without making six layers and without choosing each country manually in the Layer Styling bar panel.


texts: How to use the neighborhood modifier of grass:r.mapcalculator in PyQGIS

texts:
I would like use r.mapcalculator in a PyQGIS processing script to change values from one raster in another accordingly. I don't know how to call the variables row and col, which are now read as simple strings in the following code:

processing.runalg('grass:r.mapcalculator',{"amap": rlayer, "bmap": emptylayer,
"formula": "amap[row,col]=bmap[(%i+row*%i),(%i+col*%i)]"% (offset,scale,offset,scale),
"GRASS_REGION_PARAMETER": "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), "outfile": bdry_resamp})

Does anyone know the correct syntax to use this?

It would be fine to do it manually calling the module from the processing toolbox in QGIS as well, but then I would have to specify the 'formula' as well.


texts: Creating buffers that automatically update when new feature added using QGIS

texts:
I have currently added a CSV file of earthquakes on a 7 day basis and added in the URL so that it gets updated automatically. I also have population count data for the whole of Central and South America. What would be the best way to add in buffers (e.g. 5km, 10km, 25km), so they automatically update once an earthquake appears and how would I go about extract the population data within each buffer, so I know the population affected by each.
Moreover, I have events for the world but I only want them for Central and South America. Is there any where I can filter by attributes, so I only have events for these two areas?


texts: How to prevent partial labels from appearing in QGIS Layout Manager

texts:
How do I prevent labels from partially appearing on the QGIS print layout?

I have tried to editing the label properties and have unchecked the "Allow truncated labels on edges of map" checkbox but still this does not prevent partial labels from displaying the layout manager.

Any suggestions? I am using 3.16.4-Hannover.
I have also unchecked the "Allow truncated labels on edges of map"


texts: I am trying to create a template for a standard output map layout including all surround elements. I am doing a temporary contract for a government department. I'll be done here in a few weeks, and they don't have an in-house GIS tech. I want to leave them with a template that requires no additional work to generate a map. Within ArcGIS (10.4) layout view, elements like legend and scale bar are dynamic with respect to layers and map scale, respectively. However, I don't know if it is possible to have a grid/graticule that dynamically adjusts the interval to be appropriate to the scale of the map. Is there a way to do this?


texts: JOSM vs QGIS/uDIG for road network editing and storage

texts:
My application needs to store a road network and provide editing functionality for that network (editing/updating/removing lanes, signs, traffic lights, speed limits info, etc.).

I find JOSM tool to be very convenient - it has a really nice user experience.

This tool works with local .osm files only.

It is really easy to export the updated .osm file into a database with osm2pgrouting tool for later routing tasks.

On the other side there is QGIS/uDIG which can work with a PostGIS database layers. So, as I understand, all my road network edits will be stored in a PostGIS-enabled database.

My question is - is there any real advantage in using QGIS/uDIG for my task?

(except of getting rid of continuous usage of osm2pgrouting after JOSM edits of the .osm file).


texts: Assign field values (double) based on values of other field (string) using field calculator

texts:
I want to assign numbers to a field, WScore (double), based on field values of Soils (string). They are in the same table of a raster. I'm not trying to add text to a double. What I want to do is this:
Where Soils = Ds, Oo, Mu, P, or Se,  WScore = 2.5.

I have little python experience so any guidance would be really helpful. Here's what I've come up with so far.

Expression:

WScore=Reclass(!SOILS!)

def Reclass(Soils):

    x=1
    if (Soils == "Ds"):
        x=2.5
    elif (Soils == "Oo"):
        x = 2.5
    elif (Soils == "Mu"):
        x = 2.5
    elif (Soils == "P"):
        x = 2.5
    elif (Soils == "Se"):
        x = 2.5

I'm getting error 000539, which I'm guessing is because I can't concatenate 'str' and 'int' objects, though this error # covers many.


texts: How to avoid QGIS line width scaling when exporting geometry to KML?

texts:
I have shapefiles loaded as linear geometries in QGIS.  I use the layer properties dialog to apply some styling to them (color + linewidth).  The goal is to output each layer as KML file.  My problem is, any geometry linewidth (pixels) I enter in QGIS ends up scaled by about 2.83 factor in KML.  For example, to get a 4.0 pixels KML linewidth, I have to enter 1.4111 as width value in QGIS.

What is this related to? Misuse of my own or software bug?


texts: Unable to connect POSTGIS Table to CARTO

texts:
I am a new Carto Builder (enterprise) User and I am trying to connect a Postgres9.6.4 / PostgiS 2.3.1. table to my datasets.

Here is my script that should do the trick but I receive an extra long answer with html embeded code in it answering ERROR 500.

I tried to query the ?import/api to see if there was anything pending or processing but nothing.

Is PostgreSQL 9.6 supported ? PostGIS 2.3.1 ?

The table is 500Mo large, points, WGS 84.


texts: Return street address or lat/lon coordinates from assessor's parcel number?

texts:
In the U.S., in any arbitrary state and county (subject to availability — see diagram and description), I want to:

start with:

state
county
APN (Assessor's Parcel Number)

and return:

a physical street address (preferred)
or latitude longitude coordinates.

Is there a free way to do this on the internet?

Maybe a big-data or open source project or API perhaps? Don't worry about all data being available. Coverage areas similar to the following diagram are expected and accepted.


texts:
I have a code to transform numpy array in TIFF image using Pillow. But now I need to georeference the tiff to later use it in a GIS software. In my case, I have the latitude and longitude of each of my pixel and I would like to have my GeoTIFF tagged in a specific projection which I have the definition in a .prj along with a shapefile.

I know how to do it with GDAL but I am looking for a way to do it with Pillow or another low dependency package in Python.


texts: How to batch export layers to KML in ArcMap?

texts: Is there a way to batch export layers to KML from ArcMap (ArcGIS Desktop) by altering the layer to KML tool to have more than one layer input?


texts: Applying Maximum Likelihood Classification to a single band (NDWI)

texts:
I have applied ML to a single band file, specifically a NDWI (McFeeters version). Comparing the results of this classification, with ML applied too a normal 5 band RapidEye Scene, the first one gave me better results (the breaking waves of the Bay of Bengal weren't classified as vegetation, as in the latter).

I'd like to know why this is happening and how does the ML equation changes when it is applied to a single band file?


texts:
Our initial need is to have a GIS tool that overlays images on a map which is not facing north.  In our research we have seen that only Mapbox and OpenLayers support this.

Along with that we want other functionality - zoom, easyto , scaleto, pitch etc.

We have seen that mapbox is an answer to all problems but its pretty pricey now (499 p/month). Does that price paid for mapbox as the mapping layer or for a commercial use with mapbox.js / mapboxgl.js as well?

If mapbox.js is free to use commercially, can we use it with OpenLayers as maps layer?

We tried to use OpenLayers so far but it lacks some functionality for us, and have some issues with image georeference (for example the image must be rotated in advance to face north)

Any advice?


texts:
Working with PostgreSQL basic view in QGIS I noticed that:

Select *
from t1
where t1.column = 'some_condition'

Views doesn't support constraints according to this topic:
https://stackoverflow.com/questions/11667508/how-to-add-primary-key-to-view

My question is: Is some way to use  1:N relation QGIS with PostreSQL views?
When I try to create this relation in QGIS - Project - properties using View + table_relation: not load correctlly.

Otherwise, when I used spatial_table (with PK constraint) + table relation: works correctlly.

Someone have any idea?


texts: Unknown coordinate reference systems in UK

texts:
Can anyone identify what the following coordinate reference systems are?

This should be in Chester, UK:

This should be in Eastbourne, UK:

This should be in Horsham:


texts: Using MOD within a GROUP BY subquery on point dataset

texts:
I have a  dataset of 200 points (in shp/f gdb/p gdb format) in ArcGIS.

I am trying to use the MOD() function within a GROUP BY subquery as a definition query on the layer to thin out the data. What I need is to group the data by the 'Name' column, and then filter each 'Name' by a specific amount or percentage. Ideally I would like to use a MOD value of 5 or 10.

The table below is an example.

ID | Time | Name
---|------|-----
01 | 0100 | FOX
02 | 0200 | FOX
03 | 0300 | FOX
04 | 0400 | FOX
05 | 0600 | FOX
-
06 | 0100 | CAT
07 | 0200 | CAT
08 | 0300 | CAT
-
09 | 0100 | DOG
10 | 0200 | DOG
11 | 0300 | DOG
12 | 0400 | DOG
13 | 0500 | DOG
14 | 0600 | DOG
15 | 0700 | DOG

The below MOD works on the entire layer, but there is a risk of it skipping data like CAT where the query only shows rows 1, 5, 10, and 15 (file gdb format).

MOD("ID", 5) = 0

The below also works, but to a lesser degree. I need to thin the data, instead of just using the MIN and MAX values (personal gdb format).

[Time] in (SELECT MAX( [Time] )
FROM data_table GROUP BY [Name])

What would be the best way to achieve this? Is there a simpler way of doing it that I have overlooked?


texts:
I have some polygons and I want to calculate the inner angle of each polygon. Now, I have the geometry of every line of polygons in the same attribute table as shown in the following image.

I want to calculate the inner angle using the following equation:

inner_angle= the bearing angle of the previous line - the bearing angle of the following line +180 (in degrees).

How can I write this in the code block in the model builder and do this for each polygon respectively?

I just need code execute that :

            ( Can the iterate feature in model builder do this loop ( for each polygon)?)

for each polygon

for i = 1: last line in the polygon

if i < the last line in the polygon

**angle = the line bearing @ i - the line bearing @ (i+1) +180**

         if angle >360

            angle= angle -360

         end

        if angle < 0

            angle = angle +360

        end
     else
     angle =   the line bearing @ i - the line bearing @ first line

         if angle >360

            angle= angle -360

         end

        if angle < 0

            angle = angle +360

        end
     end

 end

The desired result as follows:


texts:
I have 20,000 polylines in a .gdb that represent about 5,000 building footprints. Some of the building footprints have pretty substantial gaps (image below). I would like to snap the polylines closed into polygons so that I can extrude the buildings.

Can this be done in GRASS?

I am hoping that there is something in the open source community with a "Clean" or some such functionality.


texts: Is it possible to create C# based plugins for QGIS (instead of Python)?

texts: My team has a lot of C# code that we've used for ArcGIS Add-in development.  We're thinking of porting some of these tools into QGIS as plugins using QT Designer.  However rewriting them with Python on the back-end isn't such a fun idea.  Is it possible to use C# on the back end instead?


texts: What are the requirements for building a Web Mapping and analysis Site?

texts:
I want to make a site like bhuvan in which it can provide different layers of information and also tools for performing GIS analysis

I know three main requirements :

client side scripting (openLayers)
Web services (a Map server) + Tiling/Caching Services
Database (PostGIS)

I have ArcGIS desktop installed on my computer...and we can digitise many of the layers.

My questions are :

How can I put these layers I created in ArcGIS in my website i.e., how to publish these maps to any open web services.
For what purpose we need Database like postGIS. How are web services differ from databases. With out database, can we build a
webGIS?


texts: Strange black pixel on edge when generating with gdaldem -hillshade

texts:
I'm getting a strange result when I'm generating an hillshade with the gdaldem hillshade command.

Here is my command (see the compute_edges parameter):

    /Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdaldem hillshade -compute_edges -alt 90 -z 8 -s 111120 /Users/jflegendre/Documents/GDAL_PROJECTS/MP/01-DEM_MERGED/out.tiff /Users/jflegendre/Documents/GDAL_PROJECTS/MP/02-HILLSHADED_DEM/hillshaded_dem.tiff

Here is a section of the output raster:

I've also tried to clip the result with a perimeter shapefile, but there still had some black pixels.

I've tried too the nearblack gdal command without effect. The image that you can see in this post is the result of the nearblack command.
So I'm stuck actually.

Can someone help me on this?

Regards


texts:
I was looking at this fiddle http://jsfiddle.net/neWnE/1/ and trying to see get attachments to work.

Since it uses the popuptemplate I should be able to add the showAttachments: true to it to get attachments to show in the popup window correct?

var template = new PopupTemplate({
        title : "Name: {NAME}",
        fieldInfos : [{
          fieldName : "ADDRESS",
          visible : true,
          label : "Address:"
        }],
        showAttachments: true
      });

I replaced the map services in the fiddle with my own that I enabled attachments on but they aren't appearing.

Is there something else that I would need to do?


texts: Starting at 1.8, there's a standalone Quantum GIS Browser application on Windows. I am unable to  locate that app on Mac. There is a 'Browser' panel for Mac but its functionality is limited.


texts: Parametric equation with latitude and longitude

texts:
I have an arc described by three points: arc center, starting point, ending point.

In order to calculate intermediate points of the arc, I would normally use the parametrized equation.

x(t) = x(center) + radius * cos(a0 + a*t)

y(t) = y(center) + radius * sin(a0 + a*t)

being a = a1 - a0; a1 = ending angle; a0 = starting angle.

t goes from 0 (starting point) to 1 (ending point).

The problem is that my points are not described in a cartesian system (x,y), but using WSG84 (latitude and longitude).

How can I calculate intermediate points like this?

My first attempt was going from LLA to ECEF like explained here (using 0 as altitude).

Then I parametrized the arc in 3D space like explained here.

I had to calculate the three vectors and parametrized equation becomes (the same of the link but in C++ and just for an arc instead of the whole circumference:

tPoint.x = ecefCenter.x + radius * (v1.x * cos(angleStart + (angleEnd - angleStart)*t) + v2.x * sin(angleStart+ (angleEnd - angleStart)*t));
tPoint.y = ecefCenter.y + radius * (v1.y * cos(angleStart+ (angleEnd - angleStart)*t) + v2.y * sin(angleStart+ (angleEnd - angleStart)*t));
tPoint.z = ecefCenter.z + radius * (v1.z * cos(angleStart+ (angleEnd - angleStart)*t) + v2.z * sin(angleStart+ (angleEnd - angleStart)*t));

I calculate arc points with it and then come back ECEF to LL.

Results are not bad but there is some approximation in any of the steps, it is not exact.

The three original points:

9:03:34.012863:          Daniel testing, turn center latitude:
9:03:34.012926:          42.4292
9:03:34.012984:          Daniel testing, turn center longitude:
9:03:34.013016:          -6.91824
9:03:34.013062:          Daniel testing, start latitude of turn:
9:03:34.013092:          42.4422
9:03:34.013127:          Daniel testing, start longitude of turn:
9:03:34.013150:          -6.89639
9:03:34.013192:          Daniel testing, end latitude of turn:
9:03:34.013218:          42.45
9:03:34.013247:          Daniel testing, end longitude of turn:
9:03:34.013267:          -6.91908

Starting and ending point, as obtained from parametrization:

9:03:34.014212:          42.4422
9:03:34.014236:          -6.89639

9:03:34.016545:          42.4499
9:03:34.016570:          -6.91662

Here you can see the problem graphically.


texts: Is it possible to avoid certain areas using osm2po? For example I would like to avoid central London. I can't see anyway to put this as a constraint in Java.


texts:
I've downloaded the SVG icon pack from https://www.mapbox.com/maki/

Using the Import from URL has never worked for me, so how do I import these symbols into QGIS?

As far as I can tell I need an .XML defining each symbol.


texts: how to transfer GPS device coordinates into mobile android app?

texts:
i'm trying to make an android app that collects data field in high accuracy, since the phone native GNSS doesn't provide highly accurate coordinates, i want to use Coordinates provided by a real GNSS device and transfer them to my app.

any way to do that using a real GNSS stand-alone device ?


texts: Converting XY Coordinates to Projected Local Coordinate System using QGIS

texts:
I am working in QGIS 3.14, and I have created a shapefile with X and Y coordinates that has COVID testing centers in Maryland and the EPSG is 4386.  I want to do a buffer algorithm on these sites to see the locations that are at least 500 meters from a current site.  However, I can only do distance in degrees and there is a warning that says I should convert to a projected local coordinate system.
How would I do so?
When I have tried to do this in the past, the points have completely disappeared from the map.


texts: How to project an image with gdalwarp to cover a Mapbox map?

texts:
I'm trying to re-project a geographical raster source covering entire world (temperature data from GFS) to Web Mercator with this command:
gdalwarp -t_srs EPSG:3857 input.grib output.tif

However, placing it as overlay over a Mapbox map shows a shift in latitude.

Original bounds: top left 0d 7'30.00"W, 90d 7'30.00"; bottom right 359d52'30.00"E, 90d 7'30.00"S
Original image dimensions: 1440x721
Projected bounds: top left 179d59'59.12"W, 81d 6'45.00"N; bottom right 179d52'20.79"E, 81d 6'38.83"S
Projected image dimensions: 1249x1016
Inspecting the output image shows that gdalwarp returned data within ~81° latitude bound, instead of expected ~85° as defined by EPSG:3857 spec. It seems this is the cause of my original issue.
What would be the correct gdalwarp params, so that the output image covers a Mapbox map?
The original file can be accessed at https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl , select any date and time, file with f000 suffix, check level 2 m above ground, variable TMP, download.


texts:
Is it possible to add two overlays in <Fill> tag.

For example, I want to add this color, #1e9e1e, and a graphic fill. xlink:href="vegetation.png".

This is my code:

            <Fill>
              <CssParameter name="fill">#1e9e1e</CssParameter>
              <GraphicFill>
                <Graphic>
                  <ExternalGraphic>
                    <OnlineResource xlink:type="simple" xlink:href="vegetation.png"/>
                    <Format>image/png</Format>
                  </ExternalGraphic>
                  <Size>15</Size>
                </Graphic>
              </GraphicFill>
            </Fill>


texts:
I've created image mosaic that does not entirely cover its bounding box and got black holes in the places without tif file for that place. I know I can use OutputTransparentColor: #000000 to make it transparent, but tif files have white color as a background. So either I have transparent color for absentee tif files or transparent color for area just around raster.

I could set black color as a background for all tifs to make OutputTransparentColor: #000000 work, but it'll take time. Is there a way to make both these colors transparent?


texts: Accessing four extent values of layer using ArcObjects in C#?

texts:
How can I access the four extent values of a layer shown in the below image using C# in ArcObjects?

Below is the code that i am using to achieve it but "AreaOfInterest" is null.

public ISpatialReference CreateSpatialRefGCS(ESRI.ArcGIS.Geometry.esriSRGeoCSType gcsType)
        {
            ISpatialReferenceFactory spatialRefFactory = new SpatialReferenceEnvironmentClass();
            IGeographicCoordinateSystem geoCS = spatialRefFactory.CreateGeographicCoordinateSystem((int)gcsType);
            return (ISpatialReference)geoCS;
        }

        public IEnvelope GetExtent(ESRI.ArcGIS.Carto.IFeatureLayer PolygonLayer)
        {
            try
            {
                IEnvelope envelope = PolygonLayer.AreaOfInterest.Envelope;
                envelope.Project(CreateSpatialRefGCS(esriSRGeoCSType.esriSRGeoCS_WGS1984));
            }
            catch (Exception ex)
            {
                return null;
            }
            return PolygonLayer.AreaOfInterest.Envelope;
        }


texts: Consolidate polygons with large spaces between them?

texts:
I am currently working in Arc 10.1 with a grid of equally spaced polygons (created with the Grid Index Features tool).  I've then broken this grid into 4 different layers, based on the number of points (form a different layer) that fall into each polygon of the grid.

This has created four patchworks, with a good deal of white space between polygons on the same layer.

What I would like to be able to do is take all the polygons in one layer and collapse them together to eliminate the white space.

Ideally, I would also like to be able to reverse the process after I've run some analyses.

I have included some images to help illustrate what I'm talking about.  (Note I don't have enough reputation to link to the larger versions of the pictures individually, but you can view a gallery of them here:

1: I have a grid of evenly spaced polygons

2: I have separated it out into four different polygon layers

3: I am working with each layer separately

4: I want to move the existing polygons on that layer alone to eliminate the empty space

5: When the empty space is eliminated, all the polygons on the layer should make a single, solid mass

6: Reversing the process would ideally return the polygons to their original orientation


texts: Distance percentage of a point relative to a line end in ArcGIS

texts:
I have a Polyline and a Points shapefile. I need to create a new Field in the Points shapefile with the distance percentage from the coordinate-priority of the line they belong to. The coordinate priority is the lowest-rightmost end of each line feature. I'm working with ArcGIS 10.2.2 (Advanced License). Here's an example:

Any suggestions?


texts: Scripting creation of cogo polylines including curves

texts:
Is it possible to create polyline features from cogo type text files? I have some files that have survey type coordinates in them and they contain straight line segments as well as curves (defined by tangent, arc and radius).

i need to do this on hundreds of files so using the cogo toolbar in arcmap is not really an option.

i know that i can parse the text files into a table or spreadsheet of coordinates and values and if the lines were straight i could figure out how to create them.

My issue is with the curves. Is there anyway to create those using arcpy, arcobjects, or sql server spatial?


texts: SAGA Raster Calculator in QGIS Graphical Modeller not subtracting

texts:
I am using QGIS 2.8.1 to build a model. All is working fine until the raster calculator input. I understand that the raster calculator in QGIS Modeller is a SAGA tool. With this, I have tried both the following formulae:

g1-g2

a-b

I keep on getting an error message, am I typing in the formula completely wrong?

All I want to do is take away one raster layer from the other, giving me a new raster output.


texts: QGIS 2: How to create a regular vector grid at fine scale and for large area?

texts:
In QGIS, we wish to create a regular vector grid with 25*25m resolution and over an area of around 12000 km2. We tried to use:
-vector>research tools> vector grid
- and v.mkgrid function from GRASS via QGIS
but this was not successful. We tried to create grid with smaller resolution and it did work. So, it seems that the issue is due to the amount of data that this procedure generates.
if anyone has an idea about how to realize such a fine resolution grid it will be really much appreciated,
thanks,
H&M


texts:
I have downloaded some SRTM Dem files here: http://dds.cr.usgs.gov/srtm/version2_1/SRTM30//srtm30
They are from the USGS site.
I would like to export them as 3d.
I have ran through a few tutorials online, e.g.

http://www.brocku.ca/maplibrary/proc...iew/3dview.htm

http://prod.library.utoronto.ca/maplib/3D.pdf
But still don't get it.

Question:

so are the SRTM DEM files exportable to be viewed in 3d via Arcscene/ Arcmap?
how do you do it?


texts: Getting height data in XML of fixed points in The Netherlands

texts:
In Holland we use a tool offered by the government to get the height data based of NAP (Nieuw Amsterdams Peil), the tool works but it's kind of clunky and inefficient.

http://www.rijkswaterstaat.nl/apps/geoservices/napinfo3/

Anyone can login using username: rws with password: nap. It's in dutch but if you zoom in on the map you will see the locations of these "bolts" and there is some data attach to each of them. Now I am wondering if I could make a HTTP GET request to get access to this data. For example, I could make an app where the user inputs a address or coordinate and it would automatically generate a map and a list with these height bolts and there data.

There are sites that seem to supply all kinds of geo data but I am having trouble finding these bolts with "exact" height data. http://geodata.nationaalgeoregister.nl/bag/wfs?service=WFS&request=GetCapabilities&. Opengis is more of a international organization and they also provide an API to access data in JSON or XML but I'm not sure they have data about these specific height bolts in The Netherlands.
http://schemas.opengis.net/gml/3.1.1/base/gmlBase.xsd

How can I get an X amount of these bolts with height in XML format through a HTTP GET Request for a random location.

edit
It seems the following should yield the server capabilities of the CKAN API: https://geodata.nationaalgeoregister.nl/napinfo/wfs?request=GetCapabilities but it results in a error.

The following URL does yield a abstract feature which is what I am actually looking to "GET" from a specific measuring point: https://geodata.nationaalgeoregister.nl/napinfo/wfs?service=WFS&version=1.1.0&request=DescribeFeatureType&

That link tells me it is possible to somehow access specific measurement bolts and show me the height, locations in the wall, the name, date of measurement, etc. I just have no idea how to create a proper GET request to get these objects near a address or coordinate.


texts: How to extract all Way elements tagged as highways from OSM data using QGIS?

texts:
I want to extract/filter in all Way elements within a specific region tagged as 'highway' from an OSM file. I am currently thinking of using a Python library called osmread (by A Dezhin) that I stumbled upon online. It's pretty straightforward. Here's an example usage that I got from his repository.

from osmread import parse_file, Way

highway_count = 0
for entity in parse_file('foo.osm.bz2'):
if isinstance(entity, Way) and 'highway' in entity.tags:
    highway_count += 1

print "%d highways found" % highway_count

This is good enough, but I was wondering if qGIS also has a similar function that allows users to filter elements (Node, Way, Relation) given a specific set of criteria. I would prefer to use qGIS for the 'instant' visualization.

In the end, I hope to build a road network and conduct some analysis.

Thanks, all, in advance.


texts:
I already have a Sentinel-2 tile imported in QGIS and a shapefile corresponding to the area I want.
How do I cut the Sentinel 2 file?


texts: Taking attributes from state road data and adding to over lapping route data

texts:
Ok so...
I am working on a project that is looking at the regional public transportation within my area and comparing it to the overlapping road systems. I want to display the different attributes that these routes run on to give recommendations as to what kind of development should go in certain areas. So the attributes that would be in play would be speed, sidewalk access, bike lane access, as well as a couple of other things.

I have created the shapefile of the routes and I am going to compare them to the ISRN shapefile that I got from DOT. How would I go about taking these ISRN attributes and moving them over the shapefiles that I have created as well as the bike paths and sidewalks.


texts:
I was wondering if I could do the following within QGIS:

I want to update a column in my point layer with the attribute information from a specific column in a polygon layer - so I have a layer of postcodes and I want to give it a column of say School Catchment and I want to populate this column with the value of the school name column in the catchment area polygon layer.

Any help would be greatly appreciated.

Thanks!

Steve


texts:
I am trying to get a count of related records (Tableform) in order to validate Editform entry (using VBS).

Any help getting myCount?

I Know there are some syntax issues here, but is this even possible in Arcpad axf?

Seems simple using VBA but something like this:

Sub Validate()
    myCount = (Select count(*) from tblRelatedRecs where ID = Pcontrol.value Group by ID
    If myCount > 5 then
        Application.MessageBox ("Huma, you erred!")
        ThisEvent.Result = False
    End If
End Sub


texts:
I have binary raster data (deforestation hotspot or not) that I would like to plot spatial correlograms for in R, to look at spatial autocorrelation over different spatial lags. I came across this website explaining several packages that would plot spatial correlograms, but they all seem to be for continuous data (to calculate Moran's I for spatial autocorrelation).

I have calculated spatial autocorrelation for my binary data using the joincount.test function in the spdep package, but cannot find any package that allows plotting of spatial correlograms using jointcount. Does anyone know of any?


texts: Tool for adding raster values to vector z axis

texts:
In ArcGIS, there is a tool called "Interpolate Shape" which takes a vector layer and underneath a raster. Then it adds the raster grid values to the vectors z-coordinates.

Is there a tool in QGIS that does the same or sth. similar?


texts: Polygon won't label out beyond 1:150,000

texts:
I am using ArcGIS 10.1.  I have a polygon that won't label zoomed out beyond 1:150,000.  I am using the Maplex Label Engine.  When I switch the Date Frame to the Standard Label Engine the label shows.  However, I need the label offset of my feature by about an inch with a leader going to it.  There is no Reference Scale in the Data Frame Properties.  The Scale Range for the layer is set to 'Show layer at all scales'.  The Labels Scale Range is set to 'Use the same scale range as the feature layer.'  My label Method is 'Label all the features the same way.'  In my Placement Properties under the Conflict Resolution tab 'Never remove (place overlapping)' is checked.  This layer is not part of any group layer.  The labels are not paused and they do not print.  The lat/long for the feature is not 0,0 and the area is not NULL.  There are no other features around it.  When I click on 'View Unplaced Labels' from the Labeling toolbar nothing happens.  I have also tried removing the shapefile from my map and adding it back again.  I have tried adding the shapefile to a new map altogether.  No luck on either.  I have even tried creating a new shapefile and feature class inside of GDB, copy and pasting the feature and tracing the outline of the feature.  Still no label.

The only way I have been able to see labels at all is when I zoom in closer than 1:150,000.

I have ended up just creating an annotation DB for the label, but still want to know why this might be happening.  Is there a way to fix it?  And/or is there a way to prevent this from happening in the future.

Update:

This is what I am trying to accomplish at a smaller, more zoomed out, scale.

Everything is showing up the way I desired and expected at this scale.  However, when I zoom out beyond this, closer to the statewide scale, this happens.

You can see the label disappears and for some reason it is not reading as an unplaced label either.


texts: How does TimeSlider Control work in OpenLayers 2?

texts:
I would like to use the TimeSlider control of Openlayers in version 2 for WMS-Time, but there is no documentation or examples.

Could someone help me understand how this control works?


texts: Broken link in preview legend style in Geoserver

texts:
Currently in Geoserver, when I attached a style to a layer in publishing tab, the preview legend style image under "default style" shows a broken image display.

I also was unable to call it directly via URL getLegendGraphic. It only downloaded an empty WMS file.

The style is working fine when I use Geoserver preview layer.
Anybody had the solution for this problem?

The Geoserver instance version is 2.5.

Edited:
Yes, after  rechecked, the WMS indeed contains information:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1" >
    <ServiceException code="MissingParameterValue" locator="service">
        Could not determine service
    </ServiceException>
</ServiceExceptionReport>

the url: http://120.10.11.18:8080/geoserver/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=image/png&WIDTH=20&HEIGHT=20&LAYER=uptb_gis_bna:lokasi_utama


texts: Can't get RTKLIB to use an external geoid

texts:
I want to write the results of processing a Rinex file using geodetic heights, and an external geoid from Spain's geographic institute IGN.

I have downloaded the geoid files from their FTP in text format from ftp://ftp.geodesia.ign.es/geoide/ascii/ with .txt and .asc extensions, and in binary GRD format (ftp://ftp.geodesia.ign.es/geoide/GRD/EGM08_REDNAP.grd), but none of them seems to work, and I am getting the exact same results as with RNX2RTKP default's geoid (EGM96, 1ºx1º), while I would expect to see a significant difference in height.

In order to tell RTKLIB that I want to use the external geoid, I run the program as follows:

rnx2rtkp.e -p 0 -k config_file.conf rinex_file.obs brdc_orbits.obs -o output_file.sol

where config_file.conf has only de following lines:

out-height         =geodetic
file-geoidfile     =geoid_filename

and geoid_filename is the .asc, .txt or .grd file.

What am I doing wrong? Am I not providing the file in the correct format, or should I set something else in the config file?


texts:
I am trying to merge 300 dataframes (these are polygons) into a single dataframe. I know how to do it two by two using rbind. Given the high number of dataframes I am looking for a way to merge them all together. Every dataframe has the same structure. One example:

Simple feature collection with 1 feature and 3 fields geometry type:
POLYGON dimension:      XY bbox:           xmin: 18.62278 ymin:
2.720688 xmax: 18.68717 ymax: 2.773334 epsg (SRID):    4326 proj4string:    +proj=longlat +datum=WGS84 +no_defs   id
source year                       geometry 1  2 HUELS POLYGON ((18.66028 2.773334...

# Load Data - EXCEL
shp_list<-(read_excel(".repository.xlsx",sheet="repos",skip=0))

# Loop to load shapefiles
for (s in 1:length(shp_list$shapefile)){

varName<-print(paste(shp_list$shapefile[s]))

  assign(varName,st_read(print(paste("./",gsub("trmc_","",shp_list$original_map[s]),"/",sep ="")),paste(shp_list$shapefile[s])))
}


texts:
In my app, I use GeoServer to serve my features and OpenLayers as a client. My data contains dates info in this format: 'yyy-mm-dd hh-mm-ss' but when it comes from GeoServer, it appears like: 2013-01-08T06:17:37+03:00 !!

How can I change its format?


texts: Latitude and Longitudes are replaced with zeros when importing excel 2007 worksheet into personal geodatabase

texts:
I am trying to import an Excel 2007 data set into a personal geodatabase so I can convert coordinate systems.  The data set includes point locations in Decimal Degrees. When I import the Excel worksheet, the latitude and longitudes in the resulting geodatabase table have been replaced with zeros.  Why is this?  This didn't occur with a different data set I have.

Thanks

Jason


texts:
I have a table with a point geometry defined as   "geom geometry(Point,28355)"
I am trying to insert data using the following:

INSERT INTO tmpwatscint ( scgid, pipeid, geom )
     SELECT s.gid, p.pipe_id, ST_AsText(ST_SetSRID(ST_LineInterpolatePoint(p.geom,ST_LineLocatePoint (p.geom,st_startpoint(s.geom))),28355)) newgeom
       FROM "tmpWatPipe" p, "tmpWatServCon" s
       WHERE s.pipe_id = p.pipe_id ;

When I run the select statement on its own, it runs fine, however, when I try to insert the data I get an error:

  ERROR:  Geometry SRID (0) does not match column SRID (28355)

I am running POSTGIS="2.3.2 r15302"


texts:
I have a VectorLayer layer with a feature collection in GeoJSON. The GeoJSON format is as following:

{
    "type": "FeatureCollection",
    "totalFeatures": "unknown",
    "features": [
        {
            "type": "Feature",
            "id": "ES.RRTN.CP.1010001",
            "geometry": {
                "type": "MultiPolygon",
                "coordinates": [
                    [
                        [
                            [
                                570389.865,
                                4722149.567
                            ],
                            [
                                570391.945,
                                4722149.769
                            ],
                            [
                                570394.939,
                                4722150.04
                            ],
                            [
                                570395.084,
                                4722144.205
                            ]
                        ]
                    ]
                ]
            },
            "properties": {
                "validTo": {
                    "simpleContent": null
                },
                "validFrom": {
                    "simpleContent": null
                },
                "inspireId": {
                    "Identifier": {
                        "localId": "1010001",
                        "namespace": "ES.RRTN.CP"
                    }
                },
                "label": "1",
                "nationalCadastralReference": "1010001",
                "areaValue": 101.2779,
                "beginLifespanVersion": {
                    "simpleContent": null
                }
            }
        }
    ]
 }

I know the feature id=1010001, so I can identify a feature by its id, localid and nationalCadastralReference.

  Is there a way to get the Feature object from the
  VectorLayer/VectorSource/Map with one of this parameters?


texts:
Is there a way to get line segments from a linestring using points with lat/lng coordinates?

ST_Line_Substring is taking float8 start/end points. Is there a way to replace the float8 with lat/lng? I have also seen ST_pointN but I have random points that are not the same as the linestrng edges. ST_Split is useful but I need to install postgis 2.0 and that will require a configuration of my whole system which is in my plan but not yet..

Regards,
A


texts: Concatenate Fields into One Field in QGIS

texts:
I have seen some posts on this, but only for ArcMap.  I am looking to take a range of values from two fields and concatenate them into one field to make a unique value.  For example:

XMAX     YMAX       GRIDID
332500   5366000    33255366
333000   5366000    33305366

In ArcMap you can do something like (Mid([XMAX]),1,4)+(Mid([YMAX]),1,4).

I know this can be done in a PostGIS database, but for the moment it needs to be done on a shapefile and I am looking to use QGIS.

Thanks


texts: Google Satellite layer does not allow zooming in and out in QGIS

texts: Google Satellite layer does not allow zooming in or out in QGIS instead the kml layers are the moving around. please help how do i fix this?


texts: Calculating EVI from Sentinel 2 imagery gives difference in results between ArcMap & QGIS?

texts:
Though EVI coefficients have not been developed for Sentinel 2 imagery, I am using the same coefficients used for MODIS to calculate EVI for Sentinel imagery.

I know that Sentinel 2 imagery is TOA reflectance and has a scaling factor of 10000. Hence, can I use the following formula to calculate EVI

(2.5* (NIR/10000-red/10000))/(NIR/10000+(6*red/10000)- (7.5*blue/10000)+1)

I did the above calculation in Arcmap and QGIS but got completely different results-
QGIS gave me a plausible range of -0.048 to 0.3915,while ArcMap included a range of -31.68 to 45.17.
The mean and std in ArcMap are 0.1718 and 0.0887 leading me to believe that these extreme values are a few pixels only. How do I clip these extra pixels of in ArcMap?


texts: Compare two Digital Elevation Model (DEM)

texts:
I would like to know what's the best method(s) to compare two DEMs without using any GUI soft like QGIS or ArcGIS.
I'd like to do it in Python or C++, with GDAL or other similar lib.
Does calculating the RMSE (Root mean square error/deviation) is useful?
How to proceed to subtract DEM1 from DEM2? And I should find a flat result, right?
Is their other ways?
Because reading for each coordinate the elevation z and compare it with the other DEM is a bit heavy I think.

Thanks for help,
eo


texts: eliminate small polygons of a PNG with gdal_sieve

texts:
I want to eliminate polygons wich consist of less than 3550 px in this PNG . I tried using gdal_sieve but there seems to be an error.

gdal_sieve -st 3550 E8.png -of PNG E8sieve3550_4.png
ERROR 6: GDALDriver::Create() ... no create method implemented for this format.

Traceback (most recent call last):
  File "C:\PROGRA~1\QGISBR~1\bin\gdal_sieve.py", line 165, in <module>
    dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
AttributeError: 'NoneType' object has no attribute 'SetGeoTransform'

I also tried "PNG" instead of PNG and changed the order of the parameters..

Edit: I tried what user30184 suggested: leave out the -of PNG part so that the output is a .tiff picture. But the result is a bit strange.

Here ist the output of gdalinfo E8.png (which I used as the source image):

Driver: PNG/Portable Network Graphics
Files: E8.png
Size is 8250, 4090
Coordinate System is `'
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4090.0)
Upper Right ( 8250.0,    0.0)
Lower Right ( 8250.0, 4090.0)
Center      ( 4125.0, 2045.0)
Band 1 Block=8250x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=8250x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=8250x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=8250x1 Type=Byte, ColorInterp=Alpha

After eliminating the small polygons I would like to use "convert E8.png -type palette -colors 256 E8-palette-256colors.png" and than gdal_polygonize..


texts: Add bands' name and description to the Metadata when stacking using rasterio

texts:
I'm using this code, of @Loïc Dutrieux, to stack 8 Landsat bands.

with rasterio.open(final_list[0]) as src0:
   meta = src0.meta

with rasterio.open('stack.tif', 'w', **meta) as dst:
for id, layer in enumerate(final_list):
    with rasterio.open(layer) as src1:

        #the .astype(rasterio.int16) forces dtype to int16
        dst.write_band(id + 1, src1.read(1).astype(rasterio.int16))

When I open the stack.tif with gdalinfo I'm not getting bands' description in the Metadata:

...
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
...
Band 1 Block=3596x1 Type=Int16, ColorInterp=Gray
NoData Value=-9999
Band 2 Block=3596x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
...

I would like something similar to this:

...
Metadata:
Band_1=Band 1 Reflectance
Band_2=Band 2 Reflectance
Band_3=Band 3 Reflectance
Band_4=Band 4 Reflectance
Band_5=Band 5 Reflectance
Band_6=Band 7 Reflectance
Band_7=Band 6 Temperature
Band_8=Fmask
Image Structure Metadata:
...
Band 1 Block=300x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
Band 2 Block=300x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
...


texts:
I have problem with this.

if I use this

   function Add_Polygon() {
     ...
       bk_map_poly_feature.setStyle(my_style);
     ...
    }

function my_style() {
    var bk_map_polygon_style = new ol.style.Style({
    text: new ol.style.Text({
        text: map.getView().getZoom() > 20 ? 'area-1' : '',
        offsetY: 10,
    })
})
    return bk_map_polygon_style;
}

It is OK. It works. But I want to area name as parameter like this.

function Add_Polygon() {
...
bk_map_poly_feature.setStyle(my_style('area-1'));
...
}

function my_style(area_name) {
var bk_map_polygon_style = new ol.style.Style({
    text: new ol.style.Text({
        text: map.getView().getZoom() > 20 ? area_name : '',
        offsetY: 10,
    })

})
return bk_map_polygon_style;

}

It gives me an error.

  TypeError: Cannot read property 'getView' of undefined

How can I solve this?


texts:
I've only used the Simple Loader in ArcCatalog.

The documentation also mentions an Object Loader which resides in ArcMap.

Where do I find this?


texts: Import 2D GML file via OGR2OGR fails with "Did not get at least 3 values or invalid number of set of coordinates"

texts:
I have a GML file that seems to mix

<gml:posList srsDimension="2">

and

<gml:posList>

records. OGR2OGR (understandably) crashes when reading these files:

ERROR 1: Did not get at least 3 values or invalid number of set of coordinates <gml:posList>205962.585 476394.629 205957.111 476398.377 205946.834 476405.243 205947.315 476406.680 205945.326 476407.938 205945.374 476407.563 205945.320 476407.086 205944.605 476405.688 205940.987 476399.969 205938.375 476395.992 205937.806 476395.734</gml:posList>

I've seen other questions on this subject (How to debug ogr2ogr ERROR 1: Did not get at least 2 values or invalid number of set of coordinates?), suggesting to add -dim 2 to the OGR2OGR commandline, but this doesn't seem to work. OSGeo4W is up-to-date, and using GDAL 1.10 / 1.11 / 2.00 / 2.01.

C:\OSGeo4W64\bin\ogr2ogr -dim 3 -f "PostgreSQL" PG:"host=dbserver user=test password=test dbname=db" mygml.gml -dim 3 -overwrite -progress

It seems the parameter is completely ignored. All objects are 2D, so this line (with -dim 3) should fail on the first object; it doesn't. Setting this to 2 also doesn't change anything. (I added two -dim parameters as a test - also doesn't change anything).

Simply skipping objects (--skipfailures) isn't really an option since I have no idea how many objects are affected.

I have about 50GB worth of GML meaning a simple open & replace job is quite hard.

Does anyone know why OGR2OGR ignores my parameter, or can perhaps offer another solution?


texts: Upload shapefile using GeoServer REST API and C#

texts:
I'd like to upload shapefiles using GeoServer REST API and C#. So far, I can create workspaces using the REST API and C#. From GeoServer documentation (specific to cURL), the difference between creating a workspace and uploading a shapefile is in the HTTP Method (workspace creation uses POST, shapefile uploading uses PUT), ContentType (workspace creation uses application/xml, shapefile uploading uses application/zip), and finally there is the --data-binary parameter in shapefile uploading.

The code I use for creating workspace is shown below. Can anyone kindly help me refactor this code so that I can use it to upload a shapefile.

        String gUrl = "http://localhost:8081/geoserver/rest/workspaces";
        WebRequest request = WebRequest.Create(gUrl);

        request.ContentType = "text/xml";
        request.Method = "POST";
        request.Credentials = new NetworkCredential("geoServer-username", "password");

        byte[] buffer = Encoding.GetEncoding("UTF-8").GetBytes("<workspace><name>new_workspace_csharp_2</name></workspace>");
        Stream requestStream = request.GetRequestStream();
        requestStream.Write(buffer, 0, buffer.Length);
        requestStream.Close();

        WebResponse response = request.GetResponse();
        Console.Write("Response from GeoServer: " + response);

EDITED:
I've managed to make some headway, and I can now upload a shapefile using C# and the GeoServer REST API. My refactoring effort, which now works, is shown below:

    public bool UploadShapeFile(string workspace, string dsName, Uri zipUri)
    {
        String fileUri = zipUri.AbsolutePath;
        Console.Write(fileUri);

        byte[] localShapeFile = readLocalShapeFile(fileUri);

        String sUrl = "http://localhost:8081/geoserver/rest/workspaces/cs_ws1/datastores/cs_ds1/file.shp";
        WebRequest request = WebRequest.Create(sUrl);

        request.ContentType = "application/zip";
        request.Method = "PUT";
        request.Credentials = new NetworkCredential("gs-username", "password");

        // byte[] buffer = Encoding.GetEncoding("UTF-8").GetBytes(@fileUri);
        Stream requestStream = request.GetRequestStream();
        requestStream.Write(localShapeFile, 0, localShapeFile.Length);
        requestStream.Close();

        WebResponse response = request.GetResponse();
        Console.Write("Response from GeoServer: " + response);

        return false;
    }

    private byte[] readLocalShapeFile(string filePath)
    {
        byte[] buffer;
        FileStream fStream = new FileStream(filePath, FileMode.Open, FileAccess.Read);
        try {
            int length = (int)fStream.Length;
            buffer = new byte[length];
            int count;
            int sum = 0;

            // Read until Read method returns 0 - End of stream reached
            while ((count = fStream.Read(buffer, sum, length - sum)) > 0)
                sum += count;
        }
        finally {
            fStream.Close();
        }

        return buffer;
    }

However, I still believe there can be a more elegant way of doing this. Any suggestions will be highly appreciated.


texts: Set min and max zoom level to include a specific layer in GeoServer vector tiles

texts: GeoServer has the capability of creating vector tiles in formats such as GeoJSON, TopoJSON or PBF from different kinds of data stores. In a test application, there is a layer group of an area which has been published in PBF format as vector tiles. The problem is I want the map to include some layers in a particular zoom level but the vector tiles contain these layers in higher zoom level which are not required. I know in GWC home page there are options to determine start and end zoom level of seeding cache tiles, but can we drop some layers in A to B zoom levels of produced PBF files?


texts: How to include ol.Attributions in openlayers 3 custom build

texts:
I am updating my custom build of OpenLayers 3 after adding some new features to my app.  Everything works fine when I call http://openlayers.org/en/v3.19.1/build/ol.js but when I attempt to create my custom build and run the app, I receive the error Error: ol.Attribution is not a constructor in the browser.

I have tried ol.Attribution, ol.Attribution.*, ol.Attribution#*, ol.Attribution, and ol.control.Attribution in the exports array of my ol-custom.json file.  Nothing seems to be working.

What do I need to add to my ol-custom.json file to make this work ?


texts: How to use TerriaJS in a Rails Application?

texts:
I am looking to install TerriaJS as a Cesium/Leaflet replacement for the geospatial piece I have running on my web application. I have looked at Github and the wiki with the TerriaJS application and I don't see a clear path as to how to install and use TerriaJS with my application.

With leaflet it's fairly easy... add a reference to the .js file and the .css file, create a div... and bob's your uncle. However, I don't see that clear path with TerriaJS. It almost looks as if TerriaJS is meant to run as a nodejs application. I see where Nationalmap has used gulp (browserify) to turn it into a running web application, however, I need to do a lot of development. I went through the examples... I cloned the repo from github, but I don't see what exactly I am "browserifying".

Does anyone know of a tutorial or a plain-English example as to how to use TerriaJS inside of a Rails (or the like) web application?


texts: QGIS Intersection tool not allowing Save file as Shapefile

texts:
In QGIS 2.16.2  I am trying to intersect two Shapefiles to create a third Shapefile, however the Intersection tool is only allowing Save file as DBR, CSV, XLSX or ODS - not SHP.  One of the Shapefiles was just made and saved using Intersection Run as Batch Process in this project.

If I start a new project and add the Shapefiles  I can intersect and create a third Shapefile.  I've intersected Shapefiles created from intersected Shapefiles  in the same project before and not had this problem.

File1.shp  intersect w/ IntersectedFile2.shp should make IntersectedFile3.shp
   - Does not work in the same project where IntersectedFile2.shp was created.
Works in new project.

Is this a bug or am I missing something?


texts: How can I update ogr2ogr to include the PostgreSQL driver?

texts:
I'm trying to use ogr2ogr to upload a shapefile to a remote PostGres installation. When I ran this command:

celenius:~ celenius$ ogr2ogr -f PostgreSQL PG:"host=255.34.00.00 user=postgres dbname=mydb password=***" Dropbox/data/roads.shp;

I got the following error message:

Unable to find driver `PostgreSQL'.
The following drivers are available:
  -> `ESRI Shapefile'
  -> `MapInfo File'
  -> `UK .NTF'
 ...
  -> `SVG'
  -> `CouchDB'
  -> `Idrisi'
  -> `ARCGEN'
  -> `SEGUKOOA'
  -> `SEGY'

(I omitted some of the drivers as they are not relevant). Is there a way that I can update ogr2ogr to include PostgreSQL? I'm using a Mac with OS 10.7.4


texts: Is there a way to pass an initial value to gdal.RasterizeLayer in Python?

texts:
I have tried gdal.RasterizeLayer(..., options=['INIT=1']),
but it's ignored. I know that some preferences like
options=['ATTRIBUTE={}'.format(att_name)]),
or burn_values=[1] can be set, so I thought there should be something similar for initial values; not to mention there's a -init val option for the gdal_rasterize commandline tool. I also tried lower case init with no success.


texts:
I'm trying to serve a view of my table after passing a GeoJSON polygon.

Te idea is get all the points inside my polygon passed as GeoJSON.

For this, I've created new SQL view, and a query:

select p.id, p.nombre
from points p, (select 1 as id, st_setSRID(st_GeomFromGeoJSON(%polygon%), 4326) as geom) as pol
where st_within(p.geom, pol.geom)

where polygon is the parameter I want to pass.

I'm getting the error ERROR: unexpected end of data (at offset 0) when I'm trying to save the query in GeoServer.
The query is tested in PostGIS and running ok.

Any Idea what could be wrong?


texts:
I am drafting an application with Leaflet. Basically, I need to add several WMS layers and then cycle through them with a timed loop. This is the code (still a draft, I don't tested it) but I am missing a function for "showLayer(i)".

  var osmUrl = 'http://{s}.tile.osm.org/{z}/{x}/{y}.png',
          osm = new L.TileLayer(osmUrl, {maxZoom: 18}),
          mymap = new L.Map('map', {layers: [osm], center: new L.LatLng(41.69906, 12.39258), zoom: 5, maxZoom:8});

    var optionsArray =['opt1', 'opt2', ...];
    var layer = new Array(myoptionsArray.length);

    for (i = 0; i < myoptionsArray.length; i++) {
        l = L.tileLayer.wms(url, {optionsArray[i]});
        layer.push = l
        l.addTo(mymap);
        }

    function myLoop(){
      ...count on i....
      showLayer(i);
      timerID = setTimeout(myLoop, 500);
    }

    myLoop();

    function showLayer(i){
        ?????
    }


texts: Shortest Path Between Points for Airborne Travel

texts: I work for an aerial surveying firm and have a point shapefile of every location we need to fly over to survey. There are approximately 500 features in this dataset. I would like to determine the optimum flight path so I may fly over each site once. I do not need/ want to use a road network. I have ArcGIS 10.3.1 and every extension.


texts: Why do I get a no module error when I try to batch process in Arc10?

texts:
I was trying to batch reproject a set of shapefiles and I kept getting the error No OS Module

As the set was relatively small I just did it individually, however I want to be able to batch process in the future.

How do I correct this problem? Any advice or answers would be greatly appreciated.

Thanks,

Jason


texts: Recently, I had a question in an exam about the approximate number of islands in the Pacific Ocean. The options were Two thousand, ten thousand, twenty thousand and thirty thousand. The answer given by the exam board was 30,000. Is this a reliable number or is there a chance that it might be 20,000 or greater than 30,000. Where can I get a somewhat scientific source to verify this ???


texts:
I have a dataset of hdr files of MYD13A2 MODIS product including 4 bands NDVI, EVI, quality, and reliability.
I want to extract the NDVI values for a group of points; I have the geographic coordinate of each point (lon, lat) using ArcMap.
How can I do it?


texts: How to save training data set [pixel values] from earth engine during classification?

texts:
Link to Earth Engine code:https://code.earthengine.google.com/4da2d54d1ece4229b76c52a9d921ca70
Below I have pasted my code:

function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = ee.Number(2).pow(10).int();
  var cirrusBitMask = ee.Number(2).pow(11).int();
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
var composite = s2.filterDate('2018-03-01', '2018-03-31')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filterBounds(ft)
                  .map(maskS2clouds)
                  .median();

// Display the results.

Map.setCenter(87.32,22.40, 6);
//Map.addLayer(composite, {bands: ['B6', 'B4', 'B3'], min: 0, max: 0.3},'fcc');

var landcover_roi = composite.clip(ft);

Map.addLayer(landcover_roi, {bands: ['B6', 'B4', 'B3'], min: 0, max: 0.3},'clip');
Map.addLayer(ft, {}, 'From Fusion Table');

// Make a palette: a list of hex strings.
//var newfc = waterbody.merge(river).merge(fallow).merge(builtup).merge(forest).merge(agriland);
var newfc=waterbody.merge(agriland).merge(built_up);
//print(newfc);
//create training data
var bands = ['B2','B3','B4','B5','B6','B7','B8', 'B8A'];
var training = landcover_roi.select(bands).sampleRegions({
  collection: newfc,
  properties: ['class'],
  scale: 30
});
print(training);
//Train the classifier
var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: 'class',
  inputProperties: bands
});
var classified = landcover_roi.select(bands).classify(classifier);
//Display classification
Map.centerObject(newfc, 8);
var trainAccuracy = classifier.confusionMatrix();
print('Resubstitution error matrix: ', trainAccuracy);
print('Training overall accuracy: ', trainAccuracy.accuracy());

//Map.addLayer(sentinel, {bands: ['B6', 'B4', 'B3'], max: 0.3}, 'sentinel image');
Map.addLayer(classified, {min: 1, max: 3, palette: ['#0f1cbf', '#3df6ff', '#dfc723','#d87cf0','#335124','#81ff41']}, 'classification');
//Map.addLayer(newfc);
//Export the image to drive Storage.
Export.image.toDrive({
  image: classified ,
  description: 'classifiedRGB',
  scale: 20,
  maxPixels:350000000,
  region:ft
});
Export.table.toDrive({
  collection: newfc,
  description:'vectorsToDriveExample',
  fileFormat: 'KML'
});


texts: ArcMap Time enabled raster catalog greyed out raster

texts:
I am having a problem trying to make an animation from a raster catalog.

In the past, I have:

Added the raster catalogue to ArcMap
Added a new date field
Enabled time on the catalogue later
Created the animation and exported

For some reason with one particular raster catalogue, when I create the animation, I get a grey box through the entire series. Each of the individual rasters has the correct information, but for some reason when the animated frames are greyed out.

Is there something I missed?

I below is a screen grab of a random raster that shows up greyed out in the animation, but definitely has data in the Raster Viewer.


texts: How to exclude some parts of raster classes based on size of area?

texts: I have a raster file with 5 classifications. I've calculated area for all of the classes via Zonal geometry, Arcmap 10. Now, I want to exclude patches which are smaller than 100 sq. km. and re-calculate the area for each class. Could anyone advise me how to do so?


texts: Graduated Style Renderer not displaying for joined csv to shp

texts: I have joined a csv file (via the add spreadsheet layer plugin) and the field types have saved correctly (i.e. they are "real" types, see attached image), and yet the graduated style renderer option is not available. I have several numeric fields showing percentages and total numbers that I would like to display as graduated classes. Any advice? Ideas on why the graduated style option is not showing? See image of current style window. I am using QGIS version 2.18.23


texts: Creating mask be created from raster (ArcGIS Pro)

texts:
I am having a difficult time with utilizing ArcPro's (Desktop) to create a mask. Please see below for data and overall output desired:
Data:
Original Raster Layer

Created Mosaic Layer

Goal: Create a heat-map mask from the Moasic 2 raster.
Details: Image 2 is a stitched mosaic created from tiles (arcpy.SplitRaster()) generated from a larger project I have created. I am trying to highlight selected tiles (i.e. the generated mosaic) as a layer on-top of the original raster. I am thinking I have to create a mask but in attempting to crated one (ArcPro->Tools->Polygon From Raster), it fails [ERROR: Raster must be of integer type].

This is what I am hoping to create:


texts:
I tried to execute ogr2ogr from python with this shell command:

ogr2ogr.exe --config PG_LIST_ALL_TABLES YES -f "PostgreSQL" -append PG:"host=hostnamet user=username password=password dbname=dbname" nasfile.xml

The quick and dirty solution with os.system doesn't work. There is no data in the database after run this script:

import os
loadfile = path\\file.xml

command = "C:\\Program Files\\QGIS Chugiak\\bin\\ogr2ogr.exe --config  PG_LIST_ALL_TABLES YES -f \"PostgreSQL\" -append PG:\"host=hostname user=username password=password dbname=dbname\" " + loadfile

os.system(command)

What is the right syntax to execute it with subprocess.call?


texts: Snapping points from layer to line in QGIS?

texts:
I need to snap the points to the lines in another layer in order to create another line shapefile.

Is there any plugin / tool available for this?
The Snap geometries doesn't take points as a reference layer.


texts: Calculating Climate Suitable Days using netCDF4 and Python?

texts:
This is a follow up to Using netCDF4 Python climate algorithm?

I’m relatively novice in using netCDF files, so I’m wondering if someone can check my logic on my process of delineating what has been defined as “Vector Genesis Days” for epidemiology research. The criteria to be met for each day are:

18°C ≤ T ≤ 32°C
 20% ≤ Q ≤ 80%
1.5 mm ≤ R ≤ 20 mm

Where,
T is Temperature
Q is Relative Humidity
R is Rainfall

I’m attempting to iterate through each of these variables by day to get a count of the number of days that meet the specified parameters per year. I want to output this final annual data (1995-2009) into centroids representative of the current netCDF file grids(if possible, I did it is a “rough” manual manner) or rasters if easier.

Some background on the data:
NCEP Reanalysis I - downloaded at [ESRL NOAA PSD][2]
Data coverage:
Study Area – Orissa (Odisha), India
2.5° by 2.5° grid … ~69 mi (111 km)
9 x 11 Arrays:
Latitude: 12N – 28N
Longitude: 74E – 94E

Variables with units:
Temperature Maximum (tmax) – K (Kelvin)
Temperature Minimum (tmin) – K(Kelvin)
Precipitation Rate (prate) – kg/m2/s
Temperature Daily Average (temp) – K(Kelvin)
Specific Humidity (shum) – kg/kg
Pressure (pres) – Pascals

A few conversions were made to create the needed daily mm/day in rainfall from precipitation rate, as well as, conversion of specific humidity to relative humidity (that’s where temperature daily average and pressure needed to be utilized). In addition, a few conversions were used for Celsius to Kelvin, etc.

My current code in Python is:

"""
    Leah Bowyer
    Python Code: Calculate Vector Genesis Days >>> Convert to CSV
    """

import numpy as np
import netCDF4 as ncdf

prate=ncdf.Dataset('N:/NCEP Reanalysis_India/prate.sfc.gauss.1995.nc','r')
tmax=ncdf.Dataset('N:/NCEP Reanalysis_India/tmax.2m.gauss.1995.nc','r')
tmin=ncdf.Dataset('N:/NCEP Reanalysis_India/tmin.2m.gauss.1995.nc','r')
shum=ncdf.Dataset('N:/NCEP Reanalysis_India/shum.2m.gauss.1995.nc','r+')  ## 'r+' when converting humidity file
temp=ncdf.Dataset('N:/NCEP Reanalysis_India/temp.2m.gauss.1995.nc','r')
press=ncdf.Dataset('N:/NCEP Reanalysis_India/pres.2m.gauss.1995.nc','r')

es=6.112 * np.exp((17.67 * (temp.variables['air'][:] - 273.15)) / ((temp.variables['air'][:] - 273.15) + 243.5))
e=shum.variables['shum'][:] *( press.variables['pres'][:] / 100)/ (0.378 * shum.variables['shum'][:] + 0.622)
rhum=(e/es) *100
print shum.variables['shum'][:] ##check values before conversion
shum.variables['shum'][:]=rhum[:]

print shum.variables['shum'][:] ##check values after conversion

##check climate parameters
a=np.logical_and(shum.variables['shum'][:] <=80, shum.variables['shum'][:] >=20)
b=np.logical_and(prate.variables['prate'][:]<= 2.31481481e-03, prate.variables['prate'][:]>= 1.73611111e-05)
c=np.logical_and(tmax.variables['tmax'][:]<=305.15, tmax.variables['tmax'][:]>=291.15)
d=np.logical_and(tmin.variables['tmin'][:]<=305.15, tmin.variables['tmin'][:]>=291.15)

##sum days satisfy daily parameter
Q = sum(a)
R = sum(b)
sum(c)
sum(d)

T=np.logical_and(c,d)
precip_satisfied=np.logical_and(a,b)
criteria_satisfied=np.logical_and(T,precip_satisfied)

sum(criteria_satisfied) ##check logic in output

##output array into .csv file with parameter qualifications (daily value) by column
output = np.column_stack((sum(T).flatten(),Q.flatten(),R.flatten(),sum(criteria_satisfied).flatten()))
np.savetxt('N:\\Python Results\\2011_India.csv',output,delimiter=',')

print sum(criteria_satisfied)

Some things to note:

I’d like to output the individual parameter results to do some controls on the overall “days” count. I’ve also only coded to handle one year at a time, but this data CAN BE downloaded continuously. Also, I’ve exported the resulting array as a .csv file that I joined to an already created centroid file. I know that this is not the most efficient way.

My code does run with somewhat reasonable results, but with some analysis it does not check out correctly.


texts: Spectral mixing model error in ImageCollection using Google Earth Engine

texts:
I would like to create an ImageCollection, in specific months for the L5 with a filter that returns only one image per year, which can be the image with less cloud coverage, then apply a spectral linear mixing model and export all processed images to the google drive.
The problem occurs when I apply the function of the spectral linear mixture model to the distinct year collection I created. I believe the problem lies in this collection. I have no experience with GEE, could someone help me?

This is what I have so far:

https://code.earthengine.google.com/b56794a03e7e4963ff226747c89518eb

var bacia = ee.FeatureCollection("users/geodenilson/Bacia_Carlos");
var roi = bacia.geometry();

// define custom function for exporting images to drive
var ExportCol = function(col, folder, scale, type,
                     nimg, maxPixels, region) {
type = type || "float";
nimg = nimg || 500;
scale = 30 || 1000;
maxPixels = maxPixels || 1e10;

var colList = col.toList(nimg);
var n = colList.size().getInfo();

for (var i = 0; i < n; i++) {
  var img = ee.Image(colList.get(i));
  var id = img.id().getInfo();
  region = roi || img.geometry().bounds().getInfo()["coordinates"];

  var imgtype = {"float":img.toFloat(),
                 "byte":img.toByte(),
                 "int":img.toInt(),
                 "double":img.toDouble()
                }

  Export.image.toDrive({
    image:imgtype[type],
    description: 'LANDSAT_5_'+id,
    folder: folder,
    //fileNamePrefix: id,
    region: region,
    scale: scale,
    maxPixels: maxPixels})
}}

var getFractionssn = function(ModeloMistura){
    // Define Endmembers
    var endmemberssn =[
            [150.0, 430.0, 270.0, 4300.0, 1500.0, 550.0], /*gv*/
            [1200.0, 2400.0, 3250.0, 4600.0, 5500.0, 3600.0], /*soil*/
            [160.0, 300.0, 170.0, 150.0, 50.0, 30.0] /*water*/
      ];
    // Uminxing Data
    var unmixedsn = (ModeloMistura)
                  .select('B1','B2','B3','B4','B5', 'B7')
                  .unmix(endmemberssn)
                  .max(0)
                  .multiply(100)
                  .byte();

  var detexsn = unmixedsn.expression(
    '120*(Soil/GV)+20', {
      'Soil': unmixedsn.select (1),
      'GV': unmixedsn.select (0)});

    return detexsn;
};

var lsCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
  .filterBounds(ee.Geometry.Point(-44.85,-18.74))
  .filter(ee.Filter.dayOfYear(182, 243))
  // Add the observation year as a property to each image.
  .map(function(img) {
    return img.set('year', ee.Image(img).date().get('year'));
  });

// Make a distinct year collection; one image representative per year.
var distinctYears = lsCol.distinct('year').sort('year');

var result = distinctYears.map(getFractionssn)

ExportCol(result, 'GEE_output', 30)


texts: How do I customize QGIS' Add Feature Widget?

texts:
I have a plugin that validates some of the data a user provides when adding a new feature by the build in Add Feature UI.

I want to disable some of the input fields or provide some default values to some of the fields, or at least let the user know which fields are required.  Is there a way to customize the Add Feature input widget?

EDIT

Maybe I wasn't clear enough in my question. I'm looking for a pyqgis solution, that is, a way to dynamically customize the input form of the Add Feature dialogue.


texts: Drawing transects using point layer as start point

texts:
I'm very new to QGIS 3 and I have a problem. I got a point.shp indicating coastal transects long the world's coast set 500 m apart. I selected a subset projecting only the sandy transects along the Dutch coastline. From these points I know want to create a 5 km long transect inland 90° to the existing shoreline. The first picture shows the available points, second picture is a zoom-in and the red lines indicate what I want to achieve.
Is there any way on how to create these inland transects perpendicular to the shore? Any suggestions, insights, tips, etc?


texts: scratch workspace in-line variable, up one level?

texts:
So I'm building a model to be published as a geoprocessing service, and I'm trying to use the 'Save to Layer File' tool to create a symobolized result of my raster output.

I've set both current and scratch workspaces to separate file gdb's, and am using in-line variables (%workspace%, %scratchworkspace%) to direct my outputs. I want to use the %scratchworkspace% in-line variable for my 'Save to Layer File' output, but to my knowledge, layer files cannot be stored in a geodatabase... therefore, I need to know how I can set the 'Save to Layer File' output to the folder that the scratch gdb is stored in.

Is it possible to specify my .lyr output to be saved "up one level" from the %scratchworkspace%??


texts: Checking whether feature class has selection using ArcPy?

texts:
In a sentence, how do you test in Arcpy whether a feature class has any selected records or not? More info below to explain:

I am doing a batch process in which a sampling layer (10km grids) is used to select polygon records (of which there are about 100,000 in total). My script returns the list of polygon IDs that intersect each 10km grid square. I then do some further processing/summarizing using that list of polygon IDs.

In instances where there are polygons within the sample square, ArcGIS performs a selection on the polygons, and it's relatively easy to iterate through those records using a searchcursor. Note I'm using v10.0.

But in instances where there are no selected records (ie there are no polygons intersecting my 10km sample square), the searchcursor would return ALL records in the database. So what's the most robust way to test programatically whether a feature class has a selection or not?

I know you can use arcpy.GetCount_management to return a number that tells you how many records are selected, and if that number = the number of records in the entire layer then you can infer that either none, or all are selected.

But that's not the same as determining if there is a selection or not. What's the best way to do that in Arcpy?


texts: What I'm doing wrong with the command r.out.ascii in GRASS GIS?

texts:
I have raster in .bil format which I want to save in GRASS GIS with r.out.ascii. I get this result in .txt file:

north: 1N
south: 0
east: 1E
west: 0
rows: 1
cols: 1
55537

I suppose it is something wrong, because of dimensions of this .txt file. I didn't change any options in r.out.ascii window.

Thank you for your help.


texts: Leaflet GeoJSON with projected coordinates (3857)

texts:
I am trying to get a GeoJSON with projected coordinates to show in Leaflet.

I have no problems loading GeoJSON with lat/lng (EPSG:4326). Here is a fiddle where I load the GeoJSON with unprojected lat/lng.

I am using the Leaflet 1.0.0-beta.2 (954e83d) in the fiddle but I am happy to use another version.

I am aware of the coordsToLatLng method of L.geoJson but no idea how to make this work.

This is the GeoJSON with the projected (EPSG:3857) coordinates.

{
    "type": "FeatureCollection",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::3857"
        }
    },
    "features": [{
        "type": "Feature",
        "properties": {

        },
        "geometry": {
            "type": "Polygon",
            "coordinates": [[
            [-625570.6, 6465993.0],
            [-305006.9, 6696510.8],
            [-546330.2, 6768547.6],
            [-445478.6, 7053093.0],
            [-279794.0, 7056694.8],
            [-359034.5, 7337638.4],
            [-211359.0, 7485313.8],
            [62380.8,   6955843.3],
            [220861.7,  6909019.4],
            [91195.5,   6700112.6],
            [213658.1,  6682103.4],
            [-625570.6, 6465993.0]]]
        }
    }]
}

Notes: I am new to Leaflet. I come from a GIS background. I have a very good understanding of projections and transformations. I have a very good understanding of OpenLayers3. I use OL3 to do things such as this on a daily basis. The data in this post is obviously just an example. I am trying to learn Leaflet in order to compare the performance with OL3. I know how to use QGIS and can unproject the data manually but this is not the point of this post.

This is not a repost. I have checked every bit of information on coordsToLatLng but didn't find any good documentation and nothing like a working example. Certainly no fiddle.

Just saying this as I have read a lot of Leaflet posts on gis.stackexchange with half-a-code. And a lot of answers which suggest to change the data. I know that Leaflet is not a GIS system and the makers are not from a GIS background. In the last hour I have just read a lot of rubbish about Leaflet using EPSG:3857 as default SRS and WGS84 being an SRS.


texts: Loading a WMS layer using standalone PyQGIS

texts:
I am using the following code to check if a WMS raster is valid.

urlWithParams = "url=crs=EPSG:4326&format=image/png&layers=tasmania&styles&url=https://demo.geo-solutions.it/geoserver/ows"
rlayer = QgsRasterLayer(urlWithParams, 'some layer name', 'wms')
rlayer.isValid()

and it return True.

But when I use the same code in my standalone pyQGIS Application. It returns False.

I don't know why it happen. Could anyone help me?
Thank very much.

Here is my whole code:

#coding=utf-8
import sys, os
import processing
from processing.core.Processing import Processing

from osgeo import gdal, osr

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5 import *

# add below environ to access qgis / qt library  #
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'C:\Program Files\QGIS 3.12\apps\Qt5\plugins'
os.environ['PATH'] += r';C:\Program Files\QGIS 3.12\apps\qgis\bin;C:\Program Files\QGIS 3.12\apps\Qt5\bin'

os.environ['GDAL_DATA'] = "C:\OSGeo4W64\share\gdal"
os.environ['GDAL_DRIVER'] = "C:\OSGeo4W64\bin\gdalplugins"
os.environ['PROJ_LIB'] = "C:\Program Files\QGIS 3.12\share\proj"
os.environ['PYTHONPATH'] = "C:\OSGeo4W64\share\gdal"
sys.path.extend([r'C:\Program Files\QGIS 3.12\apps\qgis\python',r'C:\Program Files\QGIS 3.12\apps\Python37\Lib\site-packages',r'C:\OSGeo4W64\apps\Python37\lib\site-packages\pyproj\proj_dir\share\proj'])

QgsApplication.setPrefixPath(os.environ.get("QGIS_PREFIX_PATH", "C:\Program Files\QGIS 3.12\apps\qgis"), True)
qgs = QgsApplication([], True)
qgs.initQgis()

urlWithParams = "url=crs=EPSG:4326&format=image/png&layers=tasmania&styles=&url=https://demo.geo-solutions.it/geoserver/ows"
rlayer = QgsRasterLayer(urlWithParams, 'wms layer', 'wms')
if not rlayer.isValid():
    print("Layer failed to load!")

canvas=QgsMapCanvas()
canvas.setExtent(rlayer.extent())
canvas.setLayers([rlayer])

canvas.show()

and if I use just like:

self.rlayer = QgsRasterLayer("D:\myJpg.jpg", "myJpg")

It could work.
So I don't know where is the problem.Could anyone help me?

Here is my whole code:

#coding=utf-8
import sys, os
import processing
from processing.core.Processing import Processing

from osgeo import gdal, osr

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5 import *

# add below environ to access qgis / qt library  #
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'C:\Program Files\QGIS 3.12\apps\Qt5\plugins'
os.environ['PATH'] += r';C:\Program Files\QGIS 3.12\apps\qgis\bin;C:\Program Files\QGIS 3.12\apps\Qt5\bin'

os.environ['GDAL_DATA'] = "C:\OSGeo4W64\share\gdal"
os.environ['GDAL_DRIVER'] = "C:\OSGeo4W64\bin\gdalplugins"
os.environ['PROJ_LIB'] = "C:\Program Files\QGIS 3.12\share\proj"
os.environ['PYTHONPATH'] = "C:\OSGeo4W64\share\gdal"
sys.path.extend([r'C:\Program Files\QGIS 3.12\apps\qgis\python',r'C:\Program Files\QGIS 3.12\apps\Python37\Lib\site-packages',r'C:\OSGeo4W64\apps\Python37\lib\site-packages\pyproj\proj_dir\share\proj'])

QgsApplication.setPrefixPath(os.environ.get("QGIS_PREFIX_PATH", "C:\Program Files\QGIS 3.12\apps\qgis"), True)
qgs = QgsApplication([], True)
qgs.initQgis()

urlWithParams = "url=crs=EPSG:4326&format=image/png&layers=tasmania&styles=&url=https://demo.geo-solutions.it/geoserver/ows"
rlayer = QgsRasterLayer(urlWithParams, 'wms layer', 'wms')
if not rlayer.isValid():
    print("Layer failed to load!")

canvas=QgsMapCanvas()
canvas.setExtent(rlayer.extent())
canvas.setLayers([rlayer])

canvas.show()

and if I use just like:

self.rlayer = QgsRasterLayer("D:\myJpg.jpg", "myJpg")

It could work.
So I don't know where is the problem.Could anyone help me?


texts: Creating layer with geom type wkbGeometryCollection with gdal python?

texts:
I am having trouble creating a layer with geom type as geometry collection. Thus I am unable to output to a shapefile. I have attached code below. dstLayer is null. This does not happen if geom type is polygon or multilinestring or any other.

spatialReference = osr.SpatialReference()
spatialReference.SetWellKnownGeogCS('WGS84')
driver = ogr.GetDriverByName("ESRI Shapefile")
dstPath = os.path.join("common-border", "border.shp")
dstFile = driver.CreateDataSource(dstPath)
dstLayer = dstFile.CreateLayer("layer", spatialReference, ogr.wkbGeometryCollection)


texts: ArcGIS Collector using WFS as source layer?

texts:
We have situation with external company layer service that is included in Web Portal map.
This layer service provide some environment data that changes frequently (e.g. 5 per day).
We want to have the data included in ArcGIS Collector offline mode - as a read-only, supporting layer.

But we do NOT want to copy & synch. data from external service to our own ArcGIS database several times a day, add GlobalID&Archiving (by following SDI principles we try to use referencing instead of unnecessary copying). Of course, we cannot change the external service layer (to add GlobalID etc)

Is there no way to get such layers on Collector ?


texts:
The python version of QgsPoint takes just two floating point parameters, unlike the C++ equivalent. The QGIS API docs for QgsPoint gives this example:

pt_z = QgsPoint(120, 343, 77)
pt.exportToWkt() # PointZ(120 343 77)

Here is an excerpt from the error message in a python script

point = QgsPoint( wp.longitude, wp.latitude, 0.0 )
TypeError: arguments did not match any overloaded call:
QgsPoint(): too many arguments
QgsPoint(QgsPoint): argument 1 has unexpected type 'float'
QgsPoint(float, float): too many arguments

The call works fine without the third parameter.

How do I add 3D points from a python script?


texts: Loading projection along with raster in R?

texts:
Why the operation myraster = raster("test.asc") does not load the associated projection file test.prj located in the same directory as test.asc?

Does raster() need a specific argument? It seems impractical that I have to manually set the projection when it's already defined.


texts: GeoServer how to permit to user to see only the his data

texts:
I am beginner in GeoServer use.
I want to create different users in my server, however I don't know how to allow to each user the access for his data only, and how can he grant access to other specific users. How can I do this?


texts: QgsExpression cannot get expected result when there is "\" in the string

texts:
I hope to use QgsFeatureRequest(QgsExpression) to get desired features.
The statement in QgsExpression is:
expr = QgsFeatureRequest(QgsExpression("\"planid\" in (48757) AND spi like '0\\%'")

The feature has "planid" = 48757, and "spi" field is string like '0\abcdefg' .
So I intended to match '0\\' in the QgsExpression, but it didn't work at all. I got nothing returned.
If I changed the string to "\"planid\" in (48757) AND spi like '0%'" by removing '\\', then it could return a feature.
But I really need '\\', since it's a very important character.
So how should I modify the string to make it work?
UPDATE:
From @MatthiasKuhn suggestion, I tested to use '\\\\' in QgsExpression,  like "\"planid\" in (48757) AND spi like '0\\\\%'".
It didn't work from Python scripts, but I did get a hint from it. After running the statement from "Select by Expression" in QGIS, four '\' worked fine. So this is the correct format.
Therefore, I built the statement in Python with 8 '\', and it did work.
expr = QgsFeatureRequest(QgsExpression("\"planid\" in (48757) AND spi like '0\\\\\\\\%'")
I understand why Python needs to double the number of '\', but can @MatthiasKuhn or anyone else explain why I need four '\' here in QGIS?


texts: How to create callout/balloon labels in QGIS 2.0?

texts:
I'd like to use callout/balloon labels for a vector layer in my GIS using QGIS 2.0.  The only way I can see to do this is to position an ordinary label where I want it and then draw an empty text annotation box behind the label.
I've tried using an HTML annotation but I can't work our how to link to a feature attribute using HTML - I understand that I could use [% "attribute name" %] in the label (this results in a blank label), but how do I include the layer name in the html?

I'd really appreciate some help.  My apologies if my question isn't clear - I'd be more than happy to provide further details if necessary.

Thanks


texts: Extracting roads shapefiles from old TIGER/Line RT format?

texts:
Does anyone know how to extract a roads shapefile for the US Census circa 2000?

The only file I can find from US Census (e.g. https://www2.census.gov/geo/tiger/tgrcd108/CA/tgr06037.zip for Los Angeles) is in RT format.

My questions are:

Does the RT files contain the county-based "roads" layer (all roads, described here pp3.)?
If so, which RT file is it and how to extract the roads?

Related:

https://gis.stackexchange.com/questions/21573

Also, GDAL seems to have a driver:

https://www.gdal.org/drv_tiger.html

But I can't find example code in Python or anything for using the driver to extract roads.


texts: Reading a 4 band image with rasterio (python)

texts:
I am trying to view a tif satellite image which has 4 bands. I want to remove the last band (NIR) and view the RGB image only, so I am trying to split the NIR from the rest of the image. Here is my code

import rasterio
from rasterio.plot import show
from matplotlib import pyplot as plt
from rasterio import plot
import numpy as np
#to display RGB
dataset = rasterio.open('2.tif')
%matplotlib inline
plot.show(dataset.read([1,2,3]), cmap="gray")
#to display just the red band
%matplotlib inline
plot.show(dataset.read(4), cmap="gray")

I provided a screen shot of the code and the output I am getting

Displaying just 1 band seems fine, but any idea why I keep seeing an image with a yellow and white color scheme when I try to display RGB bands together? I thought it's a cmap issue at the beginning, but even when I add 'cmap="gray"' the color of the image remains the same.


texts: Error calculating field with script exported from ModelBuilder

texts:
I am attempting to add a new field called "TimeYrCode" to a list of shape files, and then add values to this field based on the value of the "gridcode" attribute, which represents different days of the year from 1-365.

The code to calculate this field works in ModelBuilder. Exported to a python script, the code reads as follows

# ADD AND CALCULATE FIELDS TO LAYERS IN THE WORKSPACE
for fc in fcList1:
    arcpy.AddField_management(fc, "TimeYrCode", "SHORT", "4", "", "", "", "NULLABLE", "NON_REQUIRED", "")
    arcpy.CalculateField_management(fc, "TimeYrCode", "myfunc( !gridcode!)", "PYTHON", "def myfunc(gridcode):\\n      if (gridcode >= 60) and (gridcode< 152):\\n          return 1\\n      elif (gridcode >= 152) and (gridcode< 244):\\n          return 2\\n      elif (gridcode >= 244) and (gridcode< 335):\\n          return 3\\n      elif (gridcode >= 335) or (gridcode< 60):\\n          return 4\\n      else:\\n          return 9")

When I attempt to run the code in PyCharm, I obtain the following error message

def myfunc(gridcode):\n      if (gridcode >= 60) and (gridcode< 152):\n          return 1\n      elif (gridcode >= 152) and (gridcode< 244):\n          return 2\n      elif (gridcode >= 244) and (gridcode< 335):\n          return 3\n      elif (gridcode >= 335) or (gridcode< 60):\n          return 4\n      else:\n          return 9
                                                                                ^
SyntaxError: unexpected character after line continuation character

Does anyone know what could be happening here? New to python and troubleshooting these kinds of issues


texts: Error in for loop not recognizing writeLAS command

texts:
I am currently trying to run a for loop in R that creates las files from a vector of 530 objects (EABplots). I have downloaded lidR, rLiDAR, and the relevant packages that come along with them (sp, raster, etc.), but it still isn't allowing me to run my for loop. Any suggestions at all would be wonderful, thank you!

My code:

plots <- seq(1:530)
for(n in 1:length(plots))
{
  plot.n <- plots[n]
  writeLAS(EABplots[[plot.n]], filename = paste(plot.n, ".las",  sep = ""))
}

Error Message:

  Error in writeLAS(EABplots[[plot.n]], filename = paste(plot.n, ".las",  :
  unused argument (filename = paste(plot.n, ".las", sep = ""))